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The  future  Naval  Electromagnetic  Railgun  will  use  a  mega-ampere  electrical  current  to 
generate  an  electromagnetic  force  which  accelerates  a  projectile  to  hypersonic  velocities.  The 
applied  current  can  raise  the  bulk  temperature  of  the  rails  by  over  100  degrees  Celsius, 
necessitating  an  active  cooling  system  for  the  rails  to  sustain  high  rates  of  fire  without  incurring 
pennanent  damage  to  the  gun.  The  electromagnetic  force  on  the  projectile  and  the  rails  creates  a 
complicated  stress  state  that  varies  as  the  projectile  passes  along  the  rail,  first  uniaxial  then 
biaxial  compression  acts  on  the  rails. 

In  this  study,  a  system  of  cooling  channels  for  fluid  flow  down  the  length  of  the  rails  was 
considered,  and  channels  with  elliptical  cross  sections  were  examined.  Elliptical  shapes  were 
considered  due  to  the  high  surface  area  available  for  convection,  relatively  low  impact  on  the 
stress  distribution,  and  low  stress  concentration  effect.  By  treating  an  elliptical  channel  as  a 
variable  area  fin  and  varying  the  size  and  aspect  ratio  of  the  ellipse  and  the  distance  between 
channels,  the  heat  transfer  capability  of  a  channel  array  was  maximized  based  on  given  flow 
conditions  and  applied  heat  flux.  The  optimal  channel  design  was  further  constrained  by  the 
applied  compressive  stresses.  It  was  found  that  ellipses  of  different  aspect  ratios  are  optimal  for 
the  uniaxial  and  biaxial  stress  states,  and  the  optimal  channel  design  was  limited  by  the 
competing  effects  of  these  two  structural  constraints. 

In  order  to  test  the  thennal  aspect  of  the  design,  a  representative  set  of  channels  were 
machined  into  one  third  scale  copper  rails  using  wire  electrical  discharge  machining.  Tests  were 
performed  using  both  a  steady  state  heat  flux  to  detennine  the  overall  heat  transfer  coefficient 
and  transient  conditions  to  detennine  the  system  thermal  relaxation  time.  In  order  to  verify  the 
structural  aspect  of  the  design,  a  finite  element  analysis  was  done  on  the  rail  cross  section  to 
compare  the  computational  stress  concentration  factors  with  the  theoretical  conelations  used  in 
the  literature.  The  results  of  both  the  thennal  experiments  and  finite  element  analysis  were  found 
to  be  in  reasonable  agreement  with  the  predicted  results. 
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1.1  THE  NEED  FOR  A  LONG  RANGE  NAVAL  GUN 

Ever  since  the  first  cannon  was  put  on  a  ship  in  the  fifteenth  century,  the  range  of  its  guns  has 
been  a  key  factor  in  a  ship’s  warfighting  capability.  The  typical  range  of  a  shipboard  cannon  was 
1400  yards,  and  this  did  not  change  significantly  until  the  development  of  the  long  range  torpedo 
necessitated  a  longer  range  surface  weapon  in  the  early  twentieth  century  (Mel  Fisher  Maritime 
Society,  2001).  At  the  time,  numerous  advances  made  in  fire  control  systems  led  to  the 
development  of  the  early  battleship  guns  which  had  a  range  of  7000  yards  (Givens,  1999).  Long 
range  guns  provided  ships  not  only  with  the  ability  to  strike  other  ships  from  a  distance  (Sea 
Strike),  but  to  support  troops  on  shore  as  well  (Naval  Surface  Fire  Support).  The  most  effective 
long  range  gun  that  has  ever  been  used  by  the  U.S.  Navy  is  the  Mk  7  16-inch  gun  that  was 
featured  on  the  Iowa  class  battleships.  These  guns  had  a  range  of  2 1  nautical  miles  (NM)  and 
shot  projectiles  weighing  1225  kg.  However,  these  guns  and  the  Iowa  class  battleships  were 
decommissioned  for  the  last  time  after  the  first  Gulf  War.  The  Mk  45  5-inch  gun  replaced  the 
16-inch  gun  as  the  NSFS  platform  for  today’s  cruisers  and  destroyers.  The  Mk  45  has  a  range  of 
only  13  NM,  which  is  insufficient  for  the  goals  of  NSFS  and  Sea  Strike  that  have  been  set  forth 
by  the  Department  of  the  Navy  (Clark,  2002).  The  development  of  a  long  range  electromagnetic 
railgun  promises  to  fill  this  operational  need  for  the  Navy.  The  proposed  railgun  will  have  a 
range  of  more  than  200  nautical  miles,  which  not  only  meets  the  current  requirement  for  NSFS, 
but  will  redefine  the  role  of  NSFS  and  ships  in  future  conflicts.  Railguns  are  also  of  interest  to 
the  U.S.  Anny  and  U.S.  Marine  Corps  because  of  their  ability  to  achieve  the  high  velocities 
required  to  defeat  hardened  targets  such  as  tanks,  bunkers  and  buildings. 

1.2  THE  ELECTROMAGNETIC  RAILGUN 

An  electromagnetic  railgun  operates  by  sending  a  current  pulse  through  a  circuit  consisting  of 
two  conducting  rails  and  an  armature  (Figure  1).  The  current  induces  a  magnetic  field  which 
creates  a  Lorenz  force  that  accelerates  the  armature  and  the  projectile  along  the  length  of  the 
rails.  To  achieve  the  required  range  and  muzzle  velocity,  an  approximately  five  mega- ampere 
(MA)  pulse  will  be  applied  to  the  railgun  for  a  period  of  ten  milliseconds.  The  railgun  must  be 
capable  of  operating  at  a  sustained  firing  rate  of  six  to  twelve  rounds  per  minute,  which 
corresponds  to  a  pulse  every  five  to  ten  seconds  (Smith  et  al,  2005). 

Projectile 
Force  (JxB) 

Armature  \ 


Figure  1:  Schematic  showing  the  basic  principles  of  an  EM  Railgun  (Reproduced  with  permission  ofNAVSEA 

PMS  405). 
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1.3  THERMAL  MANAGEMENT  OF  A  RAIL  GUN  LAUNCHER 

The  energy  required  to  operate  a  railgun  at  a  sustained  firing  rate  can  be  stored  in  either  a 
capacitor  bank  or  in  rotating  machinery.  The  capacitor  technology  is  more  developed,  and  it  is 
currently  seen  as  the  likely  candidate  for  initial  implementation  into  the  railgun  (Bernardes  et  al, 
2003).  In  the  capacitor  system,  the  current  pulse  is  shaped  by  discharging  the  capacitors 
sequentially,  creating  as  nearly  uniform  a  pulse  as  possible.  Figure  2  shows  an  example  of  the 
anticipated  current  profile. 


Figure  2:  Current  discharge  from  a  pulse  forming  network  vs.  time.  The  colored  lines  represent  the  current  supplied 
by  individual  capacitor  banks,  and  the  black  line  represents  the  total  current.  (Reproduced  with  permission  of 

NAVSEA  PMS  405). 

The  pulse  forming  network  must  supply  close  to  150  Mega  Joules  (MJ)  of  energy  to  the  railgun 
in  a  single  current  pulse,  with  most  of  the  energy  (63  MJ)  being  transferred  to  the  kinetic  energy 
of  the  projectile  (Figure  3).  Attempts  are  being  made  to  harness  the  residual  magnetic  energy, 
while  the  residual  energy  in  the  capacitors  already  exists  in  a  usable  fonn.  However,  the  energy 
lost  due  to  resistive  heating  throughout  the  railgun  system  (30-40  MJ)  cannot  be  recovered.  This 
energy  must  be  removed  quickly  and  continuously  to  keep  the  railgun  at  a  stable  operating 
temperature.  In  a  single  pulse,  approximately  15  MJ  of  energy  is  dissipated  to  the  rails  as  heat 
(Smith  et  al,  2005).  This  has  been  found  to  increase  the  bulk  rail  temperature  by  upwards  of 
100°C.  The  localized  temperature  in  areas  that  experience  concentrated  current  can  be  much 
higher. 
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Figure  3:  Energy  distribution  from  a  single  railgun  shot.  (Reproduced  with  permission  of  NAVSEA  PMS  405). 

In  current  lab  settings,  electromagnetic  guns  are  not  fired  in  a  repetitive  mode.  Consequently, 
they  are  not  designed  with  a  thermal  management  system  (Ellis,  2004).  One  of  the  largest 
railguns  in  existence  is  located  in  Kircudbright,  Scotland.  This  railgun  is  l/3ld  the  size  of  the 
U.S.  Navy’s  notional  railgun,  and  can  only  be  fired  once  every  four  hours  because  of  the  time 
required  for  cooling  (Ellis,  2004).  In  order  for  the  railgun  to  be  fired  at  sustained  rates,  some 
type  of  active  thermal  management  system  is  necessary. 

METHODS  OF  COOLING 

The  range  of  methods  of  cooling  is  as  diverse  as  the  applications  they  are  used  on;  however,  it  is 
important  to  keep  in  mind  the  ultimate  heat  sink.  Aboard  a  ship  this  will  generally  be  sea  water, 
which  has  excellent  thermal  properties  and  a  large  heat  capacity.  Using  seawater  directly  can  be 
problematic  due  to  fouling  and  corrosion.  Therefore,  the  railgun’ s  cooling  system  is  envisioned 
as  operating  a  freshwater  loop  as  either  a  liquid  or  a  two  phase  system;  the  strengths  and 
limitations  of  these  methods  will  be  examined  below. 

LIQUID  COOLING 

Because  of  their  density,  liquids  are  a  far  more  effective  heat  transfer  medium  than  gases.  Liquid 
cooling  is  broken  into  two  major  categories,  external  and  internal  cooling. 

External  cooling  is  accomplished  by  a  fluid  transferring  heat  with  the  surface  of  an  object.  One 
effective  method  for  removing  large  amounts  of  heat  is  spray  cooling,  in  which  a  liquid  is  forced 
on  to  the  surface  of  the  hot  object  by  jet  action.  Spray  cooling  is  favorable  because  it  does  not 
involve  modifying  the  object  that  is  being  cooled.  Drawbacks  to  this  method  are  that  a  collecting 
system  for  the  runoff  water  must  be  in  place,  and  the  cooling  that  is  achieved  by  the  spray  is  not 
necessarily  uniform. 


Internal  cooling  involves  passing  a  liquid  through  some  type  of  channel  imbedded  in  the  material 
being  cooled.  Advantages  to  this  system  are  that  the  inner  part  of  an  object  can  be  directly 
cooled,  and  the  cooling  is  uniform  along  the  length  of  the  channel  as  the  liquid  is  contained  in  a 
closed  loop  process.  The  disadvantage  is  that  the  cooling  channels  reduce  the  structural  integrity 
of  the  object  being  cooled. 

TWO  PHASE  COOLING 

Excellent  heat  transfer  is  possible  from  a  system  where  a  liquid  is  allowed  to  vaporize.  The  heat 
transfer  coefficient  in  this  process  is  at  least  two  orders  of  magnitude  larger  than  the  heat  transfer 
coefficient  for  a  liquid  that  does  not  change  phase.  Difficulties  in  implementing  a  system  like 
this  include  the  high  temperatures  involved  in  vaporization  and  the  pressure  changes  that  are 
associated  with  it.  These  systems  tend  to  be  more  complex  and  difficult  to  operate  than  liquid 
cooling  systems. 

COOLING  THE  RAILGUN’S  RAILS 

Liquid  internal  cooling  has  been  chosen  by  the  Navy  as  the  most  promising  method  for  cooling 
the  rails.  Air  cooling  does  not  have  the  ability  to  transfer  heat  quickly  enough  for  railgun 
applications,  spray  cooling  was  considered  too  difficult  to  accomplish  within  an  electrically 
charged  gun  bore,  and  the  pressure  changes  associated  with  two  phase  cooling  introduce  too 
many  complications  to  the  design  process. 

1.4  FORCES  ON  THE  RAILS 

As  the  current  heats  the  rails,  it  also  causes  an  electromagnetic  force  that  acts  on  the  rails  (Figure 
4).  The  acceleration  force  on  the  projectile  can  be  calculated  from  the  current  (7)  and  the 
inductance  gradient  (L' )  of  the  rails  as  shown  below  (Bernardes  et  al,  2003).  The  same 
magnetic  pressure  that  is  accelerating  the  projectile  also  acts  to  push  the  rails  apart.  Therefore, 
the  acceleration  force  on  the  projectile  can  be  used  to  find  the  compressive  force  on  the  rails. 

F  =  —L'l2  (1) 

2 

- *  r  n  - - 
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Figure  4:  Rail  cross  section  under  uniaxial  compression. 


Once  the  armature  has  passed  a  given  point,  the  current  and  magnetic  field  transition  to  a  state 
where  they  are  distributed  around  the  rail  (Figures  5a&  5b).  The  biaxial  forces  at  this  stage  can 
be  estimated  from  the  relationship, F  =  JxB  ,  where  J is  the  linear  current  density  and  B  is  the 
magnetic  field. 
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Figure  5:  a)  Representative  magnetic  field  on  half  of  a  rail  cross  section,  b)  Linear  current  density  around  the 

perimeter  of  a  rail  (Smith  et  al,  2005) 


The  current  and  magnetic  field  are  not  uniform  on  any  given  surface.  However,  the  assumption 
of  a  uniform  biaxial  stress  state  has  been  made.  The  uniform  biaxial  stress  state  was  calculated 
based  on  the  average  current  and  magnetic  field  on  each  surface  (Figure  6). 
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Figure  6:  Uniform  biaxial  stresses  on  a  rail  cross  section. 


STRESS  CONCENTRATION  FACTORS 
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An  abrupt  change  in  the  geometry  of  an  object  can  cause  an  uneven  stress  distribution,  where  the 
stress  near  the  change  in  form  is  much  higher  than  the  stress  in  the  rest  of  the  object.  This  effect 
is  seen  near  holes,  fillets,  notches  and  even  scratches  in  a  material  (Spotts  and  Shoup,  1998).  A 
stress  concentration  factor  (K)  is  used  to  describe  the  ratio  of  local  maximum  stress  to  applied 
stress  (Eq.  2). 


Published  theoretical  equations,  experimental  measurements  and  finite  element  analysis  can  all 
be  used  to  determine  the  stress  concentration  factor  for  a  particular  geometry.  Due  to  the  ease  of 
manufacturing,  cooling  channels  typically  have  a  circular  cross  section,  which  has  a  relatively 
low  stress  concentration  factor  (K  =  3)  (Peterson,  1953). 


2.0  DESIGN  OF  AN  OPTIMAL  CHANNEL  CROSS  SECTION 
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2.1  INITIAL  ASSUMPTIONS 

For  the  purpose  of  the  optimal  channel  cross  section  design,  the  heat  load  from  the  current  will 
be  treated  as  a  steady  state  heat  flux  applied  to  the  inside  surface  of  the  rail.  This  assumption 
was  chosen  because  most  of  the  heat  developed  in  the  rails  comes  from  the  current  that  travels  on 
or  very  near  to  the  front  surface  of  the  rail  (Figure  7).  The  opposite  side  of  the  rail  is  supported 
by  an  electrically  insulating  material  and  can  be  treated  as  an  adiabatic  boundary,  where  no  heat 
transfer  occurs.  Water  was  selected  as  the  coolant  in  this  design  because  it  is  an  excellent  liquid 
for  heat  transfer  based  on  its  low  viscosity,  high  heat  capacity  and  high  thennal  conductivity. 


q"  =  0 


Figure  7:  Railgun  cross  section  with  one  insulated  side  (q"  =  0  )  and  a  constant  heat  flux  applied  to  the  other. 

Because  most  of  the  heat  is  on  one  surface  of  the  rail,  the  decision  was  made  to  align  the 
channels  in  a  single  column  close  to  the  side  where  the  heat  is  applied  (Figure  8).  The  rail  cross 
section  is  divided  into  unit  cells,  with  one  unit  cell  designated  for  each  channel  (Figure  8).  This 
subdivision  is  necessary  in  order  to  put  a  dimensional  constraint  on  the  channel  size  and  to 
ensure  the  regular  spacing  of  the  cooling  channels  within  the  cross  section. 

The  optimal  shape  for  convective  heat  transfer  must  have  a  large  surface  area  while  maintaining 
a  conduction  path  so  that  surface  area  can  be  utilized.  However,  this  may  not  coincide  with  the 
best  shape  for  maintaining  structural  integrity.  For  example,  a  rectangular  shape  has  a  large 
surface  area  but  also  a  high  stress  concentration  factor  because  of  its  comers.  In  contrast,  a  circle 
has  a  low  stress  concentration  factor  and  a  small  surface  area.  An  elliptical  cross  section  was 
chosen  for  the  cross  section  design  because  it  combines  the  high  surface  area  with  relatively  low 
stress  concentration  factors  in  the  direction  of  maximum  stress. 
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Figure  8:  Three  unit  cells  in  a  rail  cross  section  with  possible  cooling  channel  shapes. 


2.2  HEAT  TRANSFER  IN  AN  ELLIPTICAL  FIN 


A  solid  that  experiences  conduction  within  its  boundaries  and  convection  between  its  surface  and 
the  surroundings  is  called  an  extended  surface.  When  the  extended  surface  is  specifically 
designed  to  increase  the  rate  of  heat  transfer  between  the  surface  and  the  surrounding  fluid  it  is 
called  a  fin  (Incropera  and  DeWitt,  2002),  an  example  of  an  array  of  fins  is  shown  in  Figure  9. 
Fins  can  take  on  many  fonns  and  they  are  widely  used  in  applications  such  as  engine  heads,  air 
conditioners  and  computer  chips. 


Figure  9:  An  array  of  fins  used  as  a  heat  sink. 


The  thermal  perfonnance  of  a  fin  is  characterized  by  the  thermal  resistance  of  the  fin,  Rfm  which 
is  computed  as  follows: 


R 


fin 


e 

(l  fin 


(3) 


where  6  is  the  difference  between  the  temperature  of  at  the  base  of  the  fin,  Tx=o  and  the  bulk  fluid 
temperature,  7T;  and  qfm  is  the  rate  of  heat  transfer  in  the  fin. 


Because  of  symmetry,  a  half  of  a  unit  cell  can  be  analyzed  and  treated  as  a  variable  area  fin 
(Figure  10).  The  fin  with  the  lowest  Rfm  has  the  highest  rate  of  heat  transfer,  so  the  objective  of 
this  study  was  to  find  the  fin  design  (and  consequently  the  channel  design)  that  has  the  minimum 
Rfin  while  meeting  the  structural  constraints  that  will  be  later  explained. 


Figure  10:  Half  section  of  a  fin  model. 


Four  parameters  define  the  channel  shape  and  position:  the  half-length  of  the  minor  and  major 
axes  ( a  and  b,  respectively),  the  minimum  thickness  of  the  fin  along  the  y  axis  (t),  and  the 
distance  between  b  and  the  edge  of  the  unit  cell  ( d ).  These  dimensions  are  shown  in  Figure  10. 
For  the  purpose  of  analyzing  the  heat  transfer  properties  of  the  fin,  the  material  between  the 
channel  edge  and  the  surface  of  the  rail  (d)  can  be  modeled  as  a  conduction  resistance  in  series 
with  the  fin  and  does  not  need  to  be  considered  when  analyzing  the  heat  equation  on  the  interval 
-b<x<b. 


The  general  energy  equation  for  a  variable  area  fin  is  shown  as  Eq.  (4),  where  Ac  is  the  fin  cross 
sectional  area,  h  is  the  convection  heat  transfer  coefficient,  ks  is  the  thennal  conductivity  of  the 
fin,  and  As  is  the  fin  surface  area  (Incropera  and  DeWitt  2002): 
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The  energy  equation  specific  to  a  fin  with  elliptical  boundaries,  with  v*  defined  as  the  local  fin 
half  height  is  (Fisher  and  Torrance  2001): 
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y*  =  a  +  t-y  (6) 

There  is  no  heat  transfer  across  the  left  boundary  of  the  unit  cell,  so  the  boundary  condition  at 
x=-b  is  adiabatic  (Eq.  7).  The  right  boundary  condition  (x=b)  was  set  as  a  constant  temperature 
difference,  0=  Tx-Tm  (Eq.  7;  Figure  1 1). 


Figure  11:  Boundary  conditions  for  the  fin  model. 

9  is  assumed  to  be  one  and  the  rate  of  heat  transfer  (q/m)  is  then  found  for  a  unity  temperature 
difference.  Since  R  fin  is  the  ratio  of  these  two  quantities,  R  fm  is  independent  of  9  as  long  as  the 
governing  equation  is  linear  with  respect  to  9. 

qfm  is  determined  from  conduction  at  the  right  boundary  of  the  fm  (x=b)  because  at  this  point  on 
the  fin  all  of  the  thermal  resistance  comes  from  conduction. 

i  i  a  dO  ,0. 

<lfin=<lcondLb=-KAc- £;•  (g) 

This  problem  was  solved  per  unit  length  of  the  fin,  q'fin  which  reduces  to: 


f  f  I  ■  *  isl  1/ 

9fi„=qamI\x=b=-Ky 


The  boundary  value  problem  shown  in  Eqs.  (5)  -  (7)  is  a  second  order  linear  ordinary  differential 
equation,  which  was  solved  using  a  finite  difference  algorithm  in  Matlab®.  The  solutions  for  this 
equation  are  shown  in  Figure  12  for  a  channel  half  length  of  5.5  mm  ( b ),  channel  half  height  of 
2.5  mm  (a),  and  various  fin  thicknesses  ( t ). 
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Figure  12:  Nondimensional  temperature  profile  as  a  function  of  position  along  the  x  axis. 


HEAT  TRANSFER  COEFFICIENT 

An  important  variable  in  the  fin  heat  equation  is  the  convective  heat  transfer  coefficient  (h), 
which  varies  with  the  channel  geometry  and  fluid  flow  properties.  In  order  to  find  the  h  for  a 
given  channel  design,  the  geometric  parameters  a,  b  and  t  must  be  specified  as  well  as  one 
parameter  of  the  fluid  flow.  In  this  study,  a  set  pump  power  per  unit  length  for  the  total  number 
of  channels  ( W'otal )  was  used  to  detennine  the  fluid  mean  velocity.  As  the  number  of  channels 
(AO  increases,  the  power  available  for  each  channel  decreases: 


W' 

channel 


W' 

''total 
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(10) 


The  calculation  of  the  heat  transfer  coefficient  requires  the  diameter,  area  ( Aempse )  and  perimeter 
( P ellipse )  of  the  channel  (Eqs.  1 1-13). 

Ellipse  =mb  (U) 

The  perimeter  of  an  ellipse  (P ellipse)  can  not  be  defined  as  an  exact  value  without  an  infinite  series 
expansion,  so  an  approximation  was  used  (Ramanujan  1991): 

PelUpse  =  * (3(fl  +  by -  J(3a  +  b)(a  +  3b))  (12) 


Since  an  ellipse  does  not  have  a  diameter,  the  hydraulic  diameter  (D/t)  is  used  as  an 
approximation  (Eq.  13). 


(13) 
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The  Nusselt  number  (Nu)  is  a  dimensionless  parameter  that  describes  the  ratio  of  the  convective 
heat  transfer  to  the  conduction  at  the  surface  of  an  object  (Eq.  14). 


Nu  = 


hD h_ 
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(14) 


The  Nusselt  number  for  a  specific  geometry  is  a  function  of  Reynolds  number  (ReD)  and  Prandtl 
number  (Pr).  Reynolds  number  and  Prandtl  number  are  also  both  dimensionless  parameters, 
where  Reynolds  number  is  the  ratio  of  the  inertial  and  viscous  forces  acting  on  a  fluid  and 
Prandtl  number  is  the  ratio  of  kinematic  viscosity  and  thermal  diffusivity.  The  only  fluid  used  in 
this  study  was  water,  therefore  the  Prandtl  number  was  7. 152.  Several  correlations  for  the 
Nusselt  number  exist;  the  following  correlation  is  valid  for  0.5  <  Pr  <  2000  and  3000  <  ReD  < 

5(1 06),  (Incropera  and  Dewitt,  2002): 


(//8)(Reo-1000)Pr 
l  +  12.7(//8)1/2(Pr2/3-l) 


(15) 


The  pump  power  per  unit  length  ( W' )  was  set  at  a  constant  value  in  order  to  solve  for  the 
Reynolds  number  and  consequently  the  Nusselt  number  and  the  heat  transfer  coefficient. 

Because  the  pump  power  was  specified  rather  than  the  flow  rate  (V  )  or  mean  fluid  velocity  (Um), 
five  equations  had  to  be  solved  simultaneously  to  find  the  Reynolds  number.  The  following 
values  appear  in  the  equations  below:  volumetric  flow  rate  (V  ),mean  fluid  velocity  ( Um), 
pressure  loss  per  unit  length  (A p/L),  friction  factor  (/),  fluid  density  (p),  and  the  dynamic 
viscosity  ( p ).  In  Eqs.  16-20  the  following  are  unknown  values:  Reo,  Um,  V  ,  A p  ,  and / 
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2.3  STRUCTURAL  CONSTRAINTS 


As  mentioned  in  the  introduction,  both  a  uniaxial  and  a  biaxial  stress  state  are  imposed  on  the 
rails  during  each  current  pulse.  The  cooling  channel  design  must  be  able  to  withstand  loads 
under  both  conditions.  The  uniaxial  force  is  larger  than  the  biaxial  forces,  so  the  ellipse  is 
oriented  in  that  direction  to  achieve  a  low  stress  concentration  factor.  This  in  turn  makes  the 
ellipse  vulnerable  to  crushing  by  the  biaxial  stress  state. 

UNIAXIAL  FORCES 

The  uniaxial  stress  ( a  uniaxial )  imposed  on  the  rail  was  developed  from  Eq.  (2).  In  order  to  find  the 
maximum  uniaxial  stress  (Eq.  21),  the  nominal  stress  (p nominal)  was  calculated  from  the  applied 
uniaxial  stress  and  the  geometry  of  the  unit  cell  (Eq.  22),  (Figure  13). 


^"uniaxial,  max  uniaxial^] nomimal 


^  a  +  t^ 
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Figure  13:  Uniaxial  compression  on  a  unit  cell. 


For  the  uniaxial  loading  of  a  finite  width  plate  with  an  elliptical  hole,  a  third  order  polynomial 
correlation  is  used  for  the  stress  concentration  factor,  KuniwCiai ,  where  A  is  channel  width  ratio  and 
[1  is  aspect  ratio  (Young  and  Budynas,  2002). 


X  =  —  (23) 

a  +  t 

P  =  -  (24) 

a 

Knia*ial=Ci+CA  +  C3A2+C4A3  (25) 


where: 


q  =1.000 +o.oooVI7/?  +  2.000  //? 

C2  =  -0.35 1  -  0.02 17577?  -  2.483 /  j3 


C3  =3.621-5. 183VT7/?  + 4.494//? 
C4  =-2.270  +  5.204717^-4.011//? 


For  uniaxial  loading,  increasing  the  aspect  ratio  (/?)  has  the  greatest  effect  on  decreasing  the 
stress  concentration  factor  (Figure  14).  The  lowest  stress  concentration  factors  can  be  achieved 
with  a  high  aspect  ratio  and  a  low  channel  width  ratio.  For  example,  the  lowest  stress 
concentration  factor  in  Figure  14  is  found  at  the  minima  of  the  line  representing  [:>  =  3.0,  which 
corresponds  to  channel  width  ratio  of  X  =  0.3. 


Figure  14:  Uniaxial  stress  concentration  factor  vs.  X  for  various  values  of  (3. 


This  stress  concentration  factor  (Eq.  25)  is  defined  for  an  ellipse  in  a  finite  width  plate  whose 
edges  are  not  constrained  (Figure  15  a),  but  the  unit  cells  within  the  rail  actually  have  supported 
boundaries  due  to  the  symmetry  of  the  channel  (Figure  15b).  Therefore,  the  calculated  stress 
concentration  factor  for  an  object  with  free  boundaries  gives  a  conservative  estimate  of  the 
strength  of  the  unit  cell  under  uniaxial  compressive  loading. 
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a) 


Figure  15:  a)  Finite  width  plate  with  unconstrained  edges;  b)  unit  cell  with  supported  edges 


To  satisfy  the  design  constraint  due  to  the  uniaxial  stress,  the  maximum  stress  caused  by  uniaxial 
compression  ( <7uniaXiai ,  max)  must  be  less  than  the  stress  allowed  to  act  on  the  rail  ( <7aiiow )  (Eq.  26). 
The  allowable  stress  is  detennined  from  the  rail  material’s  yield  strength  and  the  Factor  of  Safety 
(. FS)  (Eq.27).  The  yield  strength  of  the  rail  material  will  affect  the  Factor  of  Safety  that  can  be 
used  in  the  design  because  a  material  with  a  higher  yield  strength  will  allow  for  a  higher  Factor 
of  Safety  while  maintaining  an  allowable  stress  that  is  greater  than  the  maximum  stress.  Because 
the  rail  material  for  the  electromagnetic  railgun  has  yet  to  be  determined,  a  baseline  Factor  of 
Safety  of  1 .0  was  used  in  this  study.  In  future  studies  the  Factor  of  Safety  can  be  raised  to  make 
the  channel  design  more  conservative  with  respect  to  structural  constraints. 
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~  ^"allow 

(26) 
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FS 

BIAXIAL  FORCES 

The  second  loading  condition  that  the  rails  are  subjected  to  is  biaxial  compression.  A  theoretical 
correlation  for  a  biaxial  stress  concentration  factor  in  a  finite  width  plane  like  the  rail  cross 
section  was  not  found  in  the  literature.  Therefore,  a  correlation  for  an  ellipse  in  an  infinite  plane 
was  examined  and  modified  to  fit  the  railgun  application. 

The  biaxial  stress  concentration  factor  for  an  ellipse  in  an  infinite  plane  varies  with  the  ratio  of 
biaxial  stresses  (a i/ai)  and  the  aspect  ratio  ([>  )  of  the  ellipse  (Figures  16  &  17). 


The  maximum  stress  occurs  at  either  point  A  or  point  B  of  the  cross  section  (Figure  16), 
depending  on  the  ratio  of  the  stresses  (Peterson  1953).  For  most  loading  conditions,  the  highest 
stresses  are  at  point  B.  For  an  ellipse  in  an  infinite  plane,  this  occurs  on  the  following  intervals: 


^<4and^>F 


CT, 
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(28) 


However,  between  the  two  intervals  given  in  Eq.  38,  the  point  of  maximum  stress  occurs  at  point 
A: 


-1  cr.  1 

P~  ^  P 


(29) 


The  maximum  stress  at  points  A  or  B  are  defined  by  the  following  equations  (Peterson,  1953): 

2 

^"biaxial  max  A  =°2(\ +-)-&<  OV  (7^  max g  =  (7,  (1  +  2/?)  -  (7,  (30) 

The  stress  concentration  factor  is  represented  as  the  ratio  of  the  maximum  stress  to  the  greater  of 
the  two  applied  stresses.  So,  depending  on  which  applied  stress  has  a  greater  magnitude,  the 
stress  concentration  factor  is  defined  as  following: 


is  _  ‘-’"biaxial,  max  is  _  ‘-’"biaxial,  max 

■“"biaxial  -“"biaxial 


cr, 


cr. 


(31) 


Due  to  the  size  of  the  cooling  channel  cross  section  relative  to  the  rail  cross  section,  the  above 
theory  (for  an  infinite  plane)  must  be  modified  in  order  to  consider  a  finite  width  plane  (Figure 
18).  This  was  achieved  by  substituting  the  nominal  stress  for  the  applied  stresses  (Eqs.  32-34). 
The  nominal  stresses  are  calculated  from  applied  stresses  on  the  rail  and  the  unit  cell  dimensions, 
as  follows: 


cr,  =  cr, 


w-2b 


2b 


a  +  t 


(32) 

(33) 
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=  4 . (1  +  2/?)-cr2 


(34) 
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Figure  18:  Biaxial  compression  acting  on  an  elliptical  channel  in  a  unit  cell. 


Since  for  all  of  the  cases  evaluated  in  this  study  oi  /  a 2  >  ///?,  the  point  of  maximum  stress 
occurred  at  point  B  (Figure  16).  The  allowable  stress  is  defined  in  the  same  way  as  for  the 
uniaxial  case  (Eq.  27,  and  the  design  constraint  for  this  stress  condition  is: 


^"biaxial  max  "^biaxial  ^*2 
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2.4  OPTIMIZATION  ROUTINE 

The  cooling  channel  design  was  optimized  for  maximum  heat  transfer  in  a  unit  cell,  which 
corresponds  to  a  minimum  thermal  resistance  (/?/„,).  This  was  accomplished  by  use  of  a 
nonlinear  constraint  minimization  routine  that  was  solved  with  respect  to  the  geometric  values  of 
the  channel  design  (a,  b,  t).  First  the  heat  equation  and  the  associated  boundary  conditions  (Eqs. 
5-7)  that  produce  Rfm  were  parameterized.  This  involved  the  equations  used  to  determine  the 
heat  transfer  coefficient  (Eq.  10-14;  15-20),  so  that  for  any  given  channel  size  the  heat  equation 
was  calculated  with  the  appropriate  heat  transfer  coefficient.  The  equations  of  maximum  stress 
acting  on  the  rail  were  also  parameterized  in  terms  of  the  channel  dimensions  (Eq.  21-25;  32-36). 
The  optimization  was  perfonned  using  a  linearly  constrained  minimization  function  in  Matlab® 
to  minimize  Rfin  while  meeting  both  the  uniaxial  and  biaxial  constraints  of  the  allowable  stress 
(Appendix  1).  Table  1  lists  all  of  the  input  variables  in  the  optimization. 


Name 

Symbol 

Units 

Cross  section  height 

Hv 

mm 

Cross  section  length 

Lx 

mm 

Number  of  channels  (unit  cells) 

N 

— 

Uniaxial  stress 

® uniaxial 

MPa 

Biaxial  Stress  1 

G  biaxial  1 

MPa 

Biaxial  Stress  2 

G biaxial  2 

MPa 

Factor  of  Safety 

FS 

1.0 

Rail  Material: 

Thermal  conductivity 

ks 

W/mK 

Yield  Strength 

°v 

MPa 

Coolant: 

Thermal  conductivity 

kf 

W/mK 

Specific  Heat  Capacity 

Cp 

J/kg-K 

Density 

P 

kg/m3 

Dynamic  viscosity 

P 

N-s/m2 

Pump  Power 

W' 

total 

W/m 

Table  1:  Input  variables  for  optimization  function 
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3.0  DESIGN  OF  AN  OPTIMAL  CHANNEL  FOR  A  FULL  SCALE  RAIL  GUN 

The  general  optimization  routine  described  in  Chapter  2  can  be  used  to  find  an  optimal  cooling 
channel  cross  section  design  for  the  full  size  rails  of  the  electromagnetic  railgun.  The  parameters 
for  the  U.S.  Navy’s  Notional  Railgun  were  used  as  a  baseline  to  for  this  design  (Table  2). 


Total  energy  supplied  to  the  Railgun 

150  Mega-Joules 

Current  needed  to  produce  this  energy 

5.3  Mega-Amperes 

Inductance  gradient  in  the  launch 

0.5  Micro-Henries  per  meter 

Force  on  the  armature 

6.25  Mega-Newtons 

Uniaxial  stress  on  the  rails 

386  Mega-Pascals 

Biaxial  stresses  on  the  rails 

240  &  180  Mega-Pascals 

Rail  cross  section 

61  Millimeters  by  135  Millimeters 

Launcher  bore  size 

135  Millimeters  by  135  Millimeters 

Table  2:  Notional  Railgun  Parameters 

3.1  INITIAL  DESIGN  CHOICES 
GEOMETRIC  CONSTRAINTS 

In  the  railgun  the  electric  current  is  concentrated  on  the  inside  surface  of  the  rail  resulting  in 
local  heating,  so  cooling  that  surface  is  the  largest  challenge  (Figure  5).  Therefore  the  decision 
was  made  to  find  a  solution  for  a  single  row  of  channels  that  are  a  specified  distance  from  the 
edge  of  the  rail  (Figure  19).  The  distance  between  the  edge  of  the  rail  and  the  tip  of  the  channel 
is  set  at  20  mm  because  most  of  the  current  is  carried  within  20  millimeters  of  the  surface  of  the 
rail  (Smith  et  al,  2005). 


20  mm 


Figure  19:  Rail  with  channels  20  mm  from  edge 


The  unit  cell  size  was  determined  based  on  the  overall  rail  cross  section  (135  mm  x  61  mm)  and 
the  number  of  channels  (AO  to  be  placed  in  the  rail  (Ellis  et  al,  2005).  Because  cooling  channels 
are  only  being  placed  along  the  height  of  the  rail  (135  mm),  as  the  number  of  cells  increased,  the 
height  ( 2a+2t )  of  the  unit  cell  decreased  while  the  width  (61mm)  remained  constant. 


MATERIAL  SELECTION 

As  calculated  from  Eqs.  (1)  and  (2),  the  initial  uniaxial  compressive  load  acting  on  the  rails  is 
386  MPa,  which  is  greater  than  the  yield  strength  of  the  standard  rail  material,  ETP  copper  (310 
MPa).  Consequently,  a  stronger  material  had  to  be  selected  just  to  design  a  rail,  regardless  of 
cooling  channel  shape.  A  beryllium  copper  alloy  (UNS  Cl 7600  TH04)  was  selected  as  the  rail 
material  for  this  study  because  of  its  high  yield  strength  (825  MPa)  and  relatively  high  thermal 
conductivity  (215  W/m-K),  (Davis  2001). 

POWER  CONSTRAINT 

The  total  pump  power  constraint  was  set  at  100  Watts  per  meter,  which  was  sufficient  to  assure 
that  there  was  always  turbulent  flow  in  the  channels.  This  power  constraint  corresponds  to 
channel  solutions  with  Reynolds  numbers  on  the  interval  23 ,000<i?ez)<7 1,000. 

3.2  RESULTS 

The  optimization  routine  was  performed  for  a  rail  cross  section  with  four  to  fourteen  unit  cells 
using  the  optimization  routine  explained  in  Chapter  2  (Appendix  2).  Figure  20  shows  that  as  the 
number  of  unit  cells  increased,  the  optimal  unit  cell  fin  resistance  increased.  It  makes  sense  that 
smaller  unit  cells  have  less  heat  transfer  ability  for  two  reasons:  (1)  the  cell  width  has  become 
thinner  relative  to  the  cell  height,  leaving  less  material  for  conduction  and  (2)  the  holes  get 
smaller  since  they  are  limited  by  the  stress  concentration  factor  they  create. 
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Figure  20:  Unit  cell  resistance  vs.  number  of  cells 

However,  the  total  resistance  of  the  rail  displays  the  opposite  trend  as  the  number  of  channels  is 
increased,  because  more  channels  provides  more  paths  for  heat  transfer.  The  rail  can  be 
represented  as  a  parallel  resistance  network.  The  total  rail  resistance  is  equal  to  the  cell 

_  ^cell 

2N  ' 


resistance  divided  by  the  number  of  fins  or  two  times  the  number  of  channels,  Rrail 
Figure  21  shows  the  total  resistance  versus  number  of  channels. 
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Figure  21:  Total  Resistance  vs.  number  of  channels. 

The  power  constraint  also  played  a  role  in  increasing  Rfm  with  an  increasing  number  of  channels. 
As  the  number  of  channels  increased,  the  amount  of  pump  power  per  channel  decreased.  This 
caused  the  mean  fluid  velocity  to  decrease.  The  decrease  in  mean  fluid  velocity  and  hydraulic 
diameter  caused  a  lower  Reynolds  number.  This  in  turn  lowered  the  Nusselt  number  which 
lowered  the  heat  transfer  coefficient.  Figure  22  shows  that  the  heat  transfer  coefficient  as  a 
function  of  the  number  of  channels. 
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Number  of  channels 

Figure  22:  Heat  transfer  coefficient  vs.  number  of  cells. 


As  the  number  of  unit  cells  increased,  the  aspect  ratio  and  the  channel  width  ratio  ( a/(t+a )) 
increased  as  well.  Figure  23a  shows  the  aspect  ratio  of  the  channel  versus  number  of  channels, 
and  Figure  23b  shows  the  channel  width  ratio  versus  number  of  channels.  In  general,  as  N 
increased,  the  optimal  ellipse  became  longer  and  had  more  area  between  each  channel  relative  to 
the  unit  cell  size.  The  aspect  ratio  of  the  channel  cross  section  changed  as  the  number  of  unit 
cells  increased  because  all  of  the  parameters  of  the  unit  cell  were  not  scaled  to  the  number  of  unit 
cells.  As  the  number  of  channels  increased,  the  cell  height  got  smaller  but  the  cell  width  stayed 
the  same.  This  caused  the  ratio  of  biaxial  stresses  to  be  different,  and  consequently  a  different 
set  of  geometric  parameters  gave  an  optimal  solution  where  both  the  uniaxial  and  biaxial 
maximum  strength  constraints  were  met.  Figure  24  shows  the  optimal  channel  design  for  four, 
eight  and  twelve  channels  in  a  rail. 
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Figure  24:  Optimal  Channel  Arrays  for  12,  8  and  4  unit  cells. 


3.3  DESIGN  CONCLUSION 

Because  of  the  magnitude  of  the  forces  imposed  on  the  rails  of  an  electromagnetic  railgun,  the 
structural  constraints  were  the  limiting  factor  in  finding  an  optimal  cooling  channel  shape. 
However,  the  heat  transfer  calculations  were  necessary  when  considering  the  number  of  channels 
to  place  in  the  rail.  Material  selection  for  the  rails  will  be  a  major  consideration  in  future  designs 
because  the  rail  material  must  have  a  high  yield  strength  as  well  as  high  electrical  and  thermal 
conductivity. 

The  presence  of  two  different  stress  states  also  drove  the  design,  because  there  are  different 
optimal  aspect  ratios  for  biaxial  and  uniaxial  stresses.  These  had  to  be  balanced  in  a  design  that 
would  meet  both  constraints. 

This  design  is  limited  by  conservative  equations  for  the  uniaxial  and  biaxial  stress  concentration 
factors.  However,  the  effect  of  neighboring  holes  on  the  stress  factor  was  not  considered.  This 
effect  is  most  severe  when  t  <  a  (which  is  not  the  case  in  this  design),  but  there  is  still  some 
interaction  when  t  >  a  (Jones  1971)  which  needs  to  be  considered. 

The  optimal  cooling  solution  for  the  rail  under  a  steady  state  heat  flux  resulted  in  the  highest 
number  of  channels  possible.  Practicality  will  limit  the  number  of  channels  in  the  rail  cross 
section,  as  well  as  the  diminishing  return  on  total  rail  resistance  as  the  channel  number  is 
increased.  Also,  the  railgun  cooling  problem  has  a  transient  nature  that  is  not  completely 
captured  by  an  applied  steady  state  heat  flux.  The  transient  behavior  needs  to  be  considered  in 
order  to  find  the  optimal  cooling  channel  design. 
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3.4  IMPLICATIONS  OF  COOLING  CHANNEL  DESIGN 

From  the  parameters  of  the  Notional  Navy  Railgun  and  the  optimized  channel  design,  the  surface 
temperature  of  the  rail  can  be  calculated  for  various  time  intervals  (Figure  25).  This  calculation 
is  based  on  both  the  thennal  resistance  of  the  fin  and  the  conduction  resistance  between  the  fin 
and  the  surface  of  the  rail  (Appendix  3). 


Figure  25:  Rail  surface  temperature  vs.  steady  state  cooling  time 

Firing  every  5.0  seconds  results  in  an  inside  surface  temperature  of  the  rails  of  -430  °C  while 
firing  every  10.0  seconds  results  in  an  inside  surface  temperature  of  ~230  °C.  These  time 
intervals  correspond  to  railgun  firing  rates  of  12  and  6  rounds  per  minute,  respectively.  The 
temperatures  at  both  of  these  firing  rates  are  high  for  sustained  operation.  However,  the  method 
used  to  calculate  the  surface  temperature  assumes  that  the  1 5  MJ  of  heat  in  the  rails  is  all  being 
applied  on  the  inside  surface  of  the  rail,  when  in  fact  the  heat  in  generated  throughout  the  rail 
cross  section.  Because  of  this  assumption,  the  surface  temperatures  predicted  here  are  higher 
than  what  would  actually  occur. 
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4.0  EXPERIMENTAL  CONFIGURATION 
4.1  SCALED  PARAMETERS 

The  proposed  rails  in  the  full  size  notional  Navy  railgun  will  be  10  to  12  meters  long  with  a  cross 
section  of  61mm  by  135  mm.  It  was  not  feasible  to  fabricate  or  heat  even  a  small  length  of  a  rail 
with  this  cross  section,  so  the  rail  size  was  scaled  down  in  order  for  a  laboratory  model  to  be 
built  and  tested.  The  experimental  rails  were  manufactured  out  of  18-inch  long  pieces  of  1-inch 
by  2-inch  Electrolytic  Tough  Pitch  (ETP)  Copper,  the  current  standard  material  for  railgun 
testing  because  of  its  high  electrical  and  thennal  conductivity  (388  W/mK,  Davis,  2001).  The 
size  reduction  in  the  rail  cross  section  corresponds  to  a  scaling  factor  of  approximately  2.5. 

Because  of  the  relatively  low  yield  strength  of  ETP  Copper  (310  MPa)  relative  to  the  uniaxial 
stress  for  full  size  rails,  the  forces  used  in  designing  the  scaled  down  rails  were  reduced  as  shown 
in  Table  3.  These  scaled  forces  are  computed  based  on  a  current  of  1  MA,  which  is  on  the  same 
order  of  magnitude  as  the  current  Army  railgun  design. 


Original  rail 

Scaled  Rail 

Ouniaxial  —  386  MPa 

r> uniaxial  110  MPa 

^biaxial  1  —  165  MPa 

Obiaxial  1  =47.1  MPa 

C>biaxial2  =  240  MPa 

^biaxial  2  —  68.6  MPa 

Table  3:  Stresses  for  a  full  size  and  scaled  rail. 


In  order  to  maintain  a  similar  Reynolds  number  with  the  design  for  the  full  size  railgun,  a 
reduced  total  pump  power  constraint  was  set  at  50.0  W/m.  This  corresponds  to  a  Reynolds 
number  of  58,600  for  the  optimal  channel  design  for  the  experimental  rails. 

In  order  for  the  test  sections  to  be  manufactured  at  a  reasonable  cost,  it  was  decided  that  only 
three  channels  should  be  placed  in  the  experimental  rail.  Manufacturing  18-inch  long  channels 
with  an  elliptical  cross  section  requires  the  use  of  Wire  Electron  Discharge  Machining  (EDM), 
which  cost  approximately  $300  for  each  channel.  In  order  for  elliptical  channels  to  be  machined 
using  this  technique,  the  wire  must  be  threaded  through  a  guide  hole.  Three  guide  holes  of  3/ 16- 
inch  diameter  were  drilled  axially  through  the  experimental  rails,  with  the  holes  evenly  spaced 
and  centered  on  the  cross  section. 


4.2  EXPERIMENTAL  CHANNEL  SHAPES 


The  optimization  code  was  run  for  the  parameters  given  above,  and  the  channel  design  shown  in 
Figure  26  was  produced. 
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0.003281  m 
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b 

0.006050  m 

Figure  26:  Optimized  Design  for  the  experimental  rail 


The  thermal  and  structural  aspects  of  this  design  were  then  independently  verified.  The 
structural  aspect  of  the  design  was  verified  using  Finite  Element  Analysis  (Chapter  5)  while  the 
heat  transfer  component  of  the  design  was  verified  by  physical  testing.  In  order  to  select  channel 
shapes  for  test,  Rfin  is  graphed  with  respect  to  b,  with  a  and  t  set  constant  for  the  values  found  in 
the  optimization  (Figure  27). 


Figure  27:  Fin  resistance  vs.  ellipse  half  length  (b). 


40 


The  optimal  solution  for  heat  transfer  occurs  at  the  minima  of  the  curve,  at  point  3.  However, 
the  optimization  code  achieves  an  increased  Rfm  when  the  design  is  limited  by  structural 
constraints  rather  than  heat  transfer  ability  (point  2).  The  channel  shapes  associated  with  both  of 
these  points  will  be  tested,  as  well  as  the  channel  design  with  a  half  length  at  point  1  (Table  4). 


Rail  1 

Rail  2 

Rail  3 

bj—  0.354  in 

b2  =  0.476  in 

b3  =  0.598  in 

Rfln  =  0.004145  K/W 

Rfln  =  0.004052  K/W 

Rfin  =  0.004028  K/W 

Non-optimized 

Channels  optimized  with 

Channels  optimized  without 

channels 

structural  constraints 

structural  constraints 

ReD=5 6,600 

ReD=58,600 

ReD=58, 700 

Table  4:  Experimental  Rail  Parameters 
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5.0  FINITE  ELEMENT  ANALYSIS  OF  DESIGN 

The  structural  aspect  of  the  channel  design  was  verified  by  comparing  stress  concentrations 
factors  from  the  finite  element  analysis  software  package  (I-DEAS)  to  stress  concentration 
factors  that  were  produced  from  the  analytical  equations.  Stress  concentration  factors  were 
compared  for  various  geometries  and  loading  conditions — from  general  cases  found  in  the 
literature  to  the  specific  rail  design  produced  by  the  optimization  code. 

For  all  of  the  I-DEAS  modeling,  a  quarter  of  a  single  unit  cell  was  used  to  represent  the  rail  cross 
section.  A  quarter  section  can  be  analyzed  as  representative  of  the  whole  rail  because  of  the 
symmetry  of  the  unit  cell  and  the  assumption  of  uniform  stresses  acting  on  the  faces  of  the  rail 
(Figure  28). 


Figure  28:  Quarter  unit  cell  with  forces  and  constraints  shown. 

The  rail  cross  sections  were  modeled  in  I-DEAS  using  a  two  dimensional  shell  mesh  of  elements 
in  plane  strain.  A  shell  mesh  in  plane  strain  assumes  that  all  defonnation  (and  consequently  all 
strain)  occurs  within  the  plane  being  modeled  (Shih,  2002).  The  assumption  that  the  part  cannot 
deform  out  of  the  plane  is  valid  for  very  thick  objects,  and  the  dimensions  of  a  rail  are  such  that 
the  rail  can  be  considered  to  have  infinite  thickness.  There  will  be  edge  effects  at  the  breach  end 
of  the  rail,  but  this  problem  will  likely  be  remedied  by  other  structural  constraints  and  this  case 
was  not  considered  in  this  design.  According  to  convention,  the  parts  were  modeled  in  the  xy 
plane  (rail  cross  section),  and  there  is  no  strain  in  the  z  direction  (along  the  length  of  the  rail). 

5.1  UNIAXIAL  STRESS  CONCENTRATION  FACTOR  FROM  I-DEAS 

The  first  subject  for  I-DEAS  modeling  was  a  general  quarter  unit  cell  with  a  circular  hole  in  it  in 
uniaxial  tension  (Figure  29). 


0.154 

0.0 
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Figure  29:  Quarter  of  a  unit  cell  with  a  circular  hole  in  uniaxial  compression,  with  stress  concentration  factors 

shown. 

The  I-DEAS  solution  displays  the  stresses  that  occur  in  each  element  in  the  model,  and  from  this 
solution  the  maximum  principal  stress  can  be  found  along  with  the  location  of  this  stress  on  the 
model.  The  uniaxial  stress  concentration  factor  was  detennined  by  finding  the  ratio  of  the  highest 
maximum  principal  stress  on  the  model  to  the  applied  stress  (Eq.  37).  For  all  uniaxial  cases,  the 
applied  stress  was  unity,  so  the  stress  concentration  factor  was  effectively  the  maximum 
principal  stress.  Also,  in  I-DEAS  compression  has  a  negative  sign,  so  the  maximum  principal 
stress  in  the  model  is  actually  the  most  negative  number  displayed  in  the  results  which 
corresponds  to  the  blue  region  on  Figure  28.  When  compared  to  the  analytical  stress 
concentration  factor,  the  difference  between  the  I-DEAS  and  analytical  stress  concentration 
factors  in  the  first  model  was  0.04%  (Eq.  37;  Appendix  4). 


K 


uniaxial 


<j. 


max  principal 


G 


applied 


(37) 


After  mesh  independence  was  verified  and  the  most  accurate  element  shape  was  found  (3  point 
triangle),  the  dimensions  of  the  plate  were  changed.  As  the  plate  got  shorter  in  the  direction  of 
the  loading  relative  to  the  hole  size,  the  numerical  results  began  to  diverge  from  the  analytical 
stress  concentration  factor.  Since  the  equation  for  the  uniaxial  stress  concentration  factor  does 
not  include  the  plate  dimension  in  the  direction  of  the  loading,  it  can  be  assumed  that  this 
equation  does  not  account  for  end  effects  in  this  direction.  The  applicability  of  the  analytical 
expression  breaks  down  under  these  conditions. 

After  the  numerical  stress  concentration  factors  for  a  circular  hole  and  their  relationship  to  the 
theory  were  well  understood,  the  stress  concentration  factor  for  an  ellipse  in  uniaixial  loading 
was  explored.  Ellipses  of  aspect  ratios  of  0.5  <  /?  <  3  were  examined  and  all  of  the  stress 
concentration  factors  were  close  to  the  analytical  results  (Appendix  4).  However,  the  stress 
concentration  factors  determined  analytically  were  consistently  higher  than  the  stress 
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concentration  factors  found  numerically.  This  is  because  the  analytical  equation  was  derived 
using  the  assumptions  of  plane  stress  rather  than  plane  strain  (Young  and  Budynas,  2002).  Plane 
stress  assumes  that  there  is  strain  but  no  stress  perpendicular  to  the  plane,  and  this  assumption  is 
valid  for  thin  plates.  Conversely,  plane  strain  models  an  object  with  infinite  thickness  under  the 
assumption  that  there  is  stress  but  no  strain  perpendicular  to  the  plane.  The  I-DEAS  modeling 
used  elements  in  plane  strain  because  these  assumptions  more  closely  match  the  physical  reality 
of  rails  in  a  railgun  even  though  the  analytical  equations  were  based  on  plane  stress.  Because 
plane  stress  gives  a  more  conservative  but  still  reasonable  stress  concentration  factor,  it  is 
acceptable  to  use  this  analytical  theory  in  the  optimization  code. 


5.2  BIAXIAL  STRESS  CONCENTRATION  FACTOR  FROM  I-DEAS 

The  correlation  between  biaxial  stress  concentration  factors  was  investigated  next.  Since  the 
analytical  solution  exists  only  for  an  infinite  plate,  assumptions  were  made  to  modify  the  theory 
for  a  finite  plate  (see  section  2.3).  Because  of  these  theoretically  sound  but  untested  assumptions 
and  the  added  variable  of  biaxial  stress  ratio,  the  biaxial  stress  concentration  factor  correlation 
proved  to  be  much  harder  to  verily  than  the  uniaxial  correlation,  which  was  not  modified  in  any 
way. 

The  first  biaxial  model  analyzed  was  a  very  large  unit  cell  with  a  very  small  hole  (hole  diameter 
1%  of  unit  cell  length  and  width)  that  closely  approximated  the  infinite  model  in  the  literature 
(Figure  30).  The  correlation  of  stress  concentration  factors  for  an  aspect  ratio  of  1.0  and  2.0 
were  almost  exact  (Figure  3 1  &  32).  The  correlation  started  to  break  down  for  an  aspect  ratio  of 
3.0  (Figure  33).  The  results  still  followed  the  trend  of  the  analytical  curve,  but  all  the  stress 
concentration  factors  were  lower  than  predicted  (Appendix  5).  This  can  be  attributed  in  part  to 
the  different  stress  theories  used  in  I-DEAS  (plane  strain)  and  the  literature  (plane  stress). 


Figure  30:  Approximation  of  an  infinite  plane  with  an  elliptical  hole  with  aspect  ratio  of  three. 


Biaxial  stress  concentration  factor,  K  .  .  , 

Diaxial 
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Figure  31:  Stress  Concentration  Factor  vs.  Stress  Ratio  for  Aspect  Ratio  of  1.0. 


Figure  32:  Stress  Concentration  Factor  vs.  Stress  Ratio  for  Aspect  Ratio  of  2.0. 
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Figure  33:  Stress  Concentration  Factor  vs.  Stress  Ratio  for  Aspect  Ratio  of  3.0. 


A  unit  cell  of  the  dimensions  that  were  found  using  the  experimental  parameters  in  the 
optimization  code  was  next  modeled  under  various  stress  ratios  (Figure  34).  The  stress 
concentration  factors  from  these  test  cases  were  compared  to  analytical  values  for  stress 
concentration  factors  that  accounted  for  the  finite  width  of  the  unit  cell  (Appendix  6).  The  stress 
concentration  factors  from  I-DEAS  followed  the  general  trend  of  the  analytical  equations,  but  the 
results  appeared  to  be  shifted  to  the  right  (Figure  35).  However,  in  the  region  of  interest  for  this 
project  (biaxial  stress  ratio  between  0.25  and  1.0),  the  correlation  between  numerical  and 
analytical  results  was  very  close,  and  the  analytical  equations  provided  a  conservative  estimate  of 
the  stress  concentration  factor. 
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Figure  34:  Finite  width  unit  cell  in  biaxial  compression  subjected  to  a  stress  ratio  of  Oi/o2=1.0,  with  stress 

concentration  factors  shown. 


Figure  35:  Stress  Concentration  factor  vs.  stress  ratio  for  an  ellipse  of  aspect  ratio  1.844  in  a  finite  width  plane. 
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5.3  FINITE  ELEMENT  ANALYSIS  CONCLUSION 

The  finite  element  analysis  perfonned  in  this  study  is  sufficient  to  verify  both  the  uniaxial  and 
biaxial  stress  concentration  factor  correlations  used  in  the  optimization  study.  The  uniaxial 
correlation  was  verified  when  there  was  sufficient  space  between  the  edge  of  the  plate  and  the 
channel,  and  the  biaxial  correlation  was  accurate  for  shapes  with  low  aspect  ratios  (J3  <  3.0)  and 
for  stress  ratios  of  -0.75<  01/02  <1.0.  All  of  these  qualifications  on  the  use  of  the  stress 
concentration  factors  fall  within  the  range  of  use  for  railgun  cooling  channels,  so  the 
optimization  routine  can  be  considered  sound  and  conservative  with  respect  to  the  structural 
constraints. 


6.0  EXPERIMENTAL  VERIFICATION  OF  HEAT  TRANSFER  PERFORMANCE 
UNDER  STEADY  STATE  CONDITIONS 
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As  explained  in  Section  3.2,  three  rails  with  different  channel  shapes  were  tested.  The  channels 
have  constant  values  of  a  and  t,  and  each  rail  has  channels  with  a  different  value  of  b  (Figure  36) 
In  Figure  35  the  four  small  holes  in  each  rail  are  taps  for  attaching  the  manifold.  Each  rail  was 
heated  and  cooled  under  steady  state  conditions,  with  its  heat  transfer  performance  measured  in 
terms  of  thermal  resistance  of  a  fin 


Figure  36:  Three  experimental  channel  designs 


6.1  EXPERIMENTAL  SET-UP 

The  heat  transfer  capability  of  the  cooling  channel  design  was  tested  by  cooling  a  rail  while  a 
constant  heat  flux  was  applied  to  its  surface  using  a  Kapton  flexible  heater  that  supplied  1000 
Watts  to  the  experimental  rail  (Figure  37). 
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Figure  37:  Kapton  heater  applied  to  rail  surface 

In  order  to  achieve  good  contact  between  the  heater  and  the  rails,  several  safeguards  were 
employed.  First,  thermal  grease  was  applied  between  the  bar  and  the  heater.  Then  the  bar  and 
heater  were  placed  in  a  containment  piece  made  of  the  composite  G-10  (Figure  38).  G-10  was 
chosen  because  it  is  fire  retardant  and  a  reasonable  insulator  (Appendix  7).  The  top  and  bottom 
piece  of  the  containment  were  bolted  together  and  a  piece  of  rubber  the  size  of  the  heater  was 
placed  under  the  heater  to  further  ensure  good  contact  between  the  heater  and  the  bar  at  all  points 
along  the  length  of  the  bar. 


Figure  38:  G-10  Containment 


Water  at  20°  C  (Toe)  was  pumped  through  the  system  using  a  Merlin  Series  Thermo  NESLAB 
Chiller.  This  chiller  was  able  to  provide  an  overall  flow  rate  ( V  )  of  approximately  1 8  liters  per 
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minute  (LPM).  Upon  exiting  the  chiller,  the  flow  was  divided  into  three  separate  lines  for  the 
three  cooling  channels.  With  the  flow  divided  and  due  to  the  losses  between  the  chiller  and  the 
flowmeters,  a  maximum  flow  rate  between  5.5  and  6.0  LPM  was  achieved  when  all  three  lines  of 
the  system  were  open.  A  flow  rate  of  6.0  LPM  corresponds  to  a  Reynolds  number  of  75,000  for 
the  optimized  channels  in  Rail  #2. 

In  order  to  pipe  the  water  from  circular  tubing  coming  from  the  chiller  into  the  elliptical 
channels,  a  manifold  device  was  designed  (Figure  39a).  Copper  tubing  was  soldered  into  one 
end  of  the  manifold,  and  then  elliptical  holes  that  were  bigger  then  all  three  experimental  channel 
sizes  were  milled  into  the  other  side  of  the  manifold.  Because  the  manifold  was  copper,  it  was 
necessary  to  insulate  it  from  the  rails.  This  was  accomplished  using  a  polyethylene  spacer  block 
with  a  rubber  gasket  on  each  side  (Figure  39b). 


51 


6.2  DATA  ACQUISITION 


The  temperature  data  was  collected  from  both  the  water  and  the  copper  rail.  Inlet  and  outlet 
water  temperatures  were  measured  by  thermocouples  inserted  into  the  flow  path,  and  the  copper 
temperature  was  measured  by  thennocouples  place  along  the  length  of  the  rail  at  one  inch 
intervals.  These  thermocouples  were  placed  to  read  the  temperature  centered  in  between  two  of 
the  channels  and  even  with  the  edge  of  the  channels  (Figure  40). 


Thermocouple  hole 


Figure  40:  a)  diagram  showing  the  point  where  the  thermocouple  measured  temperature 


The  thermocouples  used  in  the  experiment  were  calibrated  at  both  the  freezing  point  and  boiling 
point  of  water,  and  this  correction  was  applied  to  the  thennocouple  readings  (Appendix  8). 


The  thermocouples  were  read  using  a  SCXI-1000  National  Instruments  Data  Acquisition  Unit, 
with  a  SCXI-1303  Voltage  Input  Module.  The  raw  data  was  then  recorded  as  temperature  values 
using  the  software  package  Labview™.  Data  was  collected  at  flow  rates  from  2.0  to  5.5  LPM 
for  the  three  rails  with  different  values  of  b  (Table  5). 


2.0  LPM 

2.5  LPM 

3.0  LPM 

3.5  LPM 

4.0  LPM 

4.5  LPM 

5.0  LPM 

5.5  LPM 

Rail  #1  (bj) 

X 

X 

X 

X 

X 

X 

X 

X 

Rail  #2  (b2) 

X 

X 

X 

X 

X 

X 

X 

X 

Rail  #3  (b3) 

X 

X 

X 

X 

X 

X 

X 

X 

Table  5:  Test  matrix,  X  =  test  performed. 


6.3  RESULTS 


The  data  was  graphed  in  terms  of  the  difference  from  the  average  inlet  water  temperature,  Txin  , 

where  the  average  water  temperature  was  calculated  using  steady  state  temperatures  from  all 
three  channels  (Appendix  9).  Tx(z),  The  local  water  temperature  along  the  length  of  the  bar  was 
interpolated  based  on  the  inlet  and  outlet  water  temperatures  under  the  assumption  that  the  water 
temperature  increases  linearly  along  the  length  of  the  bar  (Eqs.  38-40). 


T  = 

oo  .in 


(T  -1-  T  J.T  1 

V  water, in,  1  water, in, 2  water, in, 3/ 


T 

CO, out 


(T  +T  -\-T  t 

V  water, out  A  water, out, 2  water,  out,  3  J 


(38) 

(39) 
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W  =  JV  +  z  r-*}  (40) 

As  the  cooling  fluid  enters  the  channel,  both  thermal  and  velocity  boundary  layers  begin  to 
develop  along  the  surface  of  the  channel.  These  boundary  layers  grow  down  the  length  of  the 
channel  until  they  eventually  meet  in  the  center  of  the  channel,  at  which  point  the  flow  is 
considered  to  be  fully  developed.  Before  this  point,  the  flow  is  in  the  entry  region,  and  the 
length  of  this  region  is  called  the  entry  length.  The  velocity  and  thermal  entry  lengths  are  usually 
different,  and  both  are  a  function  of  the  fluid  and  channel  properties  (Incropera  and  DeWitt, 
2002). 

The  thermally  fully  developed  region  is  important  in  analyzing  these  results  because  when  the 
flow  is  thermally  fully  developed,  both  the  temperature  difference  between  the  fluid  and  the 
surface  and  the  convective  heat  transfer  coefficient  are  constant.  Therefore,  data  points  in  this 
region  can  be  compared  to  the  predicted  values  calculated  in  the  optimization  sequence,  because 
the  heat  transfer  coefficient  used  in  the  optimization  sequence  was  based  on  fully  developed  flow 
conditions. 

In  the  fully  developed  region  the  measured  rail  temperature  should  change  at  the  same  rate  as  the 
fluid  temperature  (with  respect  to  z).  Under  this  assumption,  a  line  parallel  to  the  fluid 
temperature  was  applied  to  the  rail  temperature  data  at  the  corresponding  flow  rate,  and  the  data 
points  that  fit  the  slope  of  this  line  were  designated  as  the  fully  developed  region  (Figures  41- 
43).  These  data  points  were  chosen  by  visual  inspection,  but  because  of  the  small  number  of 
data  points  and  the  distance  between  them  it  is  likely  that  a  statistical  correlation,  such  as  a  least 
squares  fit,  would  yield  the  same  results.  In  an  ideal  experiment  the  fully  developed  region 
should  extend  to  the  end  of  the  rail,  but  this  trend  is  not  observed  due  to  heat  losses  and  changes 
in  the  flow  pattern  near  the  ends  of  the  rail. 

As  the  channel  size  increases  from  rails  #1  to  #3,  the  length  required  for  the  flow  to  become  fully 
developed  increases,  and  the  size  of  fully  developed  region  decreases  (Table  6;  Figures  41-43; 
Appendix  10).  This  is  as  predicted,  because  the  general  calculation  of  thermal  entry  length  (z ^ 
t)  is  a  function  of  the  hydraulic  diameter  of  the  channel  for  turbulent  flow  (Eq.  4 1 ;  Incropera  and 
DeWitt,  2002). 


10< 


Zfd,  T 


<60 


(41) 


Flow  rate  (LPM) 

2.0 

2.5 

3.0 

3.5 

4.0 

4.5 

5.0 

5.5 

Rail  #1 

7-13 

6-13 

6-13 

6-13 

6-13 

5  -  13 

5  -  13 

5  -  13 

Rail  #2 

9-  12 

9-13 

8-15 

8  -  15 

9-15 

9-15 

9-15 

9-15 

Rail  #3 

12-15 

10-15 

10-15 

10-15 

9-15 

8-15 

8-15 

8-15 

Table  6:  Thermally  fully  developed  regions  (in  z)  for  given  flow  rates  and  channel  geometries. 
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- 2  LPM:  water 

°  2  LPM:  copper 

- 2.5  LPM:  water 

°  2.5  LPM:  copper 
3  LPM:  water 

3  LPM:  copper 
3.5  LPM:  water 

°  3.5  LPM:  copper 

4  LPM:  water 
4  LPM:  copper 

- 4.5  LPM:  water 

°  4.5  LPM:  copper 

- 5  LPM:  water 

°  5  LPM:  copper 

- 5.5  LPM:  water 

°  5.5  LPM:  copper 


Z  (in) 


Figure  41:  Rail  #1:  Temperature  Difference  vs.  Position  along  the  length  of  the  rail 


- 2  LPM:  water 

°  2  LPM:  copper 

- 2.5  LPM:  water 

°  2.5  LPM:  copper 
3  LPM:  water 

3  LPM:  copper 
3.5  LPM:  water 

°  3.5  LPM:  copper 

4  LPM:  water 
4  LPM:  copper 

- 4.5  LPM:  water 

°  4.5  LPM:  copper 

- 5  LPM:  water 

°  5  LPM:  copper 

- 5.5  LPM:  water 

°  5.5  LPM:  copper 


Figure  42:  Rail  #2:  Temperature  Difference  vs.  Position  along  the  length  of  the  rail 
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region 


—  2  LPM:  water 
°  2  LPM:  copper 

—  2.5  LPM:  water 
°  2.5  LPM:  copper 

3  LPM:  water 

3  LPM:  copper 
3.5  LPM:  water 

0  3.5  LPM:  copper 

4  LPM:  water 
4  LPM:  copper 

—  4.5  LPM:  water 

0  4.5  LPM:  copper 

—  5  LPM:  water 

°  5  LPM:  copper 

—  5.5  LPM:  water 

0  5.5  LPM:  copper 


0  2  4  6 

Figure  43:  Rail  #3:  Temperature  Difference  vs.  Position  along  the  length  of  the  rail 


The  most  general  trend  observed  from  these  three  graphs  is  that  higher  flow  rates  are  more 
effective  in  cooling  the  rail,  and  that  an  increased  flow  rate  increases  the  heat  transfer 
performance  of  the  cooling  system  more  than  the  changes  in  channel  design. 

At  low  flow  rates,  the  smaller  channels  cool  the  rail  more  effectively  (rail  #1  was  the  best  and 
rail  #3  was  the  worst).  At  high  flow  rates,  it  is  clear  that  the  channels  in  rail  #2  and  rail  #3  are 
more  effective  that  those  in  rail  #1  but  it  is  difficult  to  determine  from  visual  inspection  if  bar  2 
or  bar  three  is  better. 

There  are  anomalies  in  the  graphs  that  need  to  be  explored.  In  rail  #1,  the  temperature  data  is 
inconsistent.  This  can  be  attributed  to  the  fact  that  the  holes  for  the  thermocouples  were  not 
drilled  to  a  precise  depth.  This  problem  was  fixed  for  the  following  experiments  and  the  data  for 
rails  #2  and  #3  show  much  more  consistent  results.  In  the  data  collected  from  rail  #2,  the  data  at 
high  flow  rates  from  the  thennocouple  at  z= 8  inches  had  to  be  removed  because  the 
thermocouple  moved  during  the  course  of  the  experiment. 

6.4  DATA  ANALYSIS 

The  local  overall  heat  transfer  coefficient  per  unit  length,  UA'(z)  was  detennined  at  each  of  the 
17  points  where  there  was  temperature  data.  UA'(z)  was  calculated  from  the  constant  amount  of 
heat  applied  to  the  rails  per  unit  length  (  Q1 )  and  the  difference  between  the  local  rail 
temperature,  Trau(z)  and  the  local  fluid  temperature  Tj/z),  where  A  T(z)  =  Traii(z)-T,,  (z)  (Eq.  42). 
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The  local  fin  resistance,  Rfm( z)  is  the  reciprocal  of  the  local  overall  heat  transfer  coefficient  (Eq. 
43).  The  mean  fin  resistance,  Rj  was  determined  by  averaging  the  individual  fin  resistances 

for  the  thermally  fully  developed  region  (Eq.  44).  The  results  are  graphed  in  Figure  44  for  each 
rail  at  each  flow  rate  (Appendices  9-10). 
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The  trends  in  Figure  44  show  that  Rail  #  1  has  the  lowest  Rfm  values  and  Rail  #  3  has  the  highest 
Rfm  values  for  a  given  flow  rate.  This  appears  to  be  the  opposite  of  what  was  predicted  in  Figure 
27  and  Table  4,  where  the  channels  in  Rail  #  3  are  expected  to  be  the  optimum  design  for  heat 
transfer.  However,  for  a  given  volumetric  flow  rate  (V)  and  increasing  values  of  channel  cross 
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sectional  area  (Aempse),  the  mean  fluid  velocity  (Um)  decreases  (Eq.  17).  To  accurately  compare 
the  performance  of  the  three  rail  designs  Rfm  must  be  analyzed  with  respect  to  Um  instead  of 
volumetric  flow  rate  (Figure  45;  Appendix  10).  For  a  given  mean  fluid  velocity  (e.g.  0.8  m/s), 
the  lowest  values  of  Rfm  correspond  to  Rail  #  3  and  the  highest  values  correspond  to  Rail  #  1 . 


As  shown  in  Figures  44  and  45,  the  Rfm  data  followed  the  predicted  Rfin  values  to  a  degree  that 
this  experiment  can  be  considered  to  validate  the  theory  used  in  the  optimization  routine.  The 
experimental  Rfin  values  follow  the  general  trend  of  the  predicted  results,  and  all  the  experimental 
values  are  within  20  %  of  the  predicted  values. 

In  the  cases  where  the  actual  Rfm  values  were  lower  than  predicted,  this  can  be  attributed  to  the 
fact  that  the  heat  does  not  actually  flow  unidirectionally  as  assumed  in  the  optimization.  Because 
heat  actually  flows  three  dimensionally,  there  are  more  paths  for  the  heat  to  dissipate  than 
assumed.  Another  reason  for  lower  Rfm  values  than  predicted  is  that  while  attempts  were  made  to 
insulate  the  rail,  no  material  is  truly  nonconductive,  so  heat  is  lost  to  the  surroundings.  When  the 
surroundings  are  not  considered  in  the  calculations,  the  heat  lost  makes  the  Rfm  values  lower  than 
they  should  be. 
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7.0  TRANSIENT  CONSIDERATIONS 

Although  the  optimization  used  in  this  study  was  based  on  a  steady  state  analysis,  the  railgun 
cooling  problem  is  inherently  transient  in  nature.  Further  optimization  and  design  will  need  to  be 
done  with  transient  assumptions,  but  as  a  starting  point,  this  steady  state  design  was  tested  to 
evaluate  its  transient  cooling  ability. 

7.1  FINITE  DIFFERENCE  MODEL  OF  THE  TRANSIENT  THERMAL  RESPONSE 

The  transient  problem  was  modeled  by  analyzing  a  half  unit  cell  broken  into  nodes  and  elements 
(Figure  46).  The  heat  transfer  was  modeled  using  an  implicit  1-D  finite  difference  solution. 


Figure  46:  Unit  cell  half  section,  with  nodes  and  elements  identified 

An  energy  balance  is  performed  on  each  element.  Figure  47a  shows  heat  flowing  into  node  i 
from  nodes  i-1  and  i+1  as  well  as  from  the  water  in  the  channel  (Too).  Because  of  the  symmetry 
of  a  unit  cell,  no  heat  crosses  the  upper  unit  cell  boundary,  so  the  top  of  the  unit  cell  can  be 
represented  as  being  insulated.  Several  geometric  parameters  need  to  be  defined  to  solve  this 
problem;  they  are  shown  in  Figure  47b. 


Figure  47:  Single  node  in  contact  with  the  cooling  channel,  with  a)  nodal  temperatures  and  heat  terms;  b)  geometric 

parameters  defined 

An  array  was  developed  for yj,  v/  and  LS;  where  LS  is  the  length  exposed  to  convection,  y}  is  the 
right-hand  boundary  of  the  element,  and  v/  in  the  average  height  of  the  element.  The  variables 
defined  above  can  be  substituted  into  the  general  heat  equation  for  a  transient  system  (Eq.  45). 


The  equation  is  then  manipulated  into  the  form  shown  in  Eq.  47,  with  the  corresponding 
coefficients  a,  b,  c  and  d. 
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When  the  heat  equation  is  generated  for  each  node  in  the  channel  and  the  corresponding  a,  b,  c 
and  d  values  are  identified,  a  tridiagonal  matrix  of  all  of  the  terms  can  be  generated  (Eq.  48). 
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The  matrix  represents  a  system  of  i  equations  with  i  unknowns  (temperature  values  at  each 
node),  so  it  can  be  solved  simultaneously  to  find  the  unknown  values  using  a  numerical 
algorithm  for  solving  a  tridiagonal  matrix  (Press  et  al,  1992).  In  order  to  find  the  temperature 
profile  of  the  bar  at  a  specific  time,  the  matrix  in  Eq.  48  must  be  generated  for  each  time  step  and 
the  values  from  the  previous  time  step  inserted  as  part  of  the  d  matrix  (Figure  48;  Appendix  1 1). 
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:  Temperature  profile  of  the  channel  cross  section  at  various  timesteps  in  a  1.0  second  interval. 


—  0.83  seconds 

450 

—  0.40  seconds 

1.67  seconds 

3.33  seconds 

6.67  seconds 

400 

— 10.0  seconds 

350 


50  100  150  200  250  300  350  400 

Node  Number  (x-direction) 


Figure  48b):  Temperature  profile  of  the  channel  cross  section  at  various  timesteps  in  a  10.0  second  interval. 
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The  average  final  temperature  of  the  bar  is  computed  using  a  weighted  average  to  account  for  the 
difference  in  volume  between  the  elements  (Eq.  49).  The  results  of  the  transient  analysis  can  be 
nondimensionalized  into  a  ratio  of  Q,  the  total  energy  transferred  from  the  bar  over  time  t,  to  Q(), 
the  maximum  possible  energy  transfer.  The  numerical  results  are  displayed  for  times  from  1  to 
10  seconds  in  Table  7. 


T 

bar,  final 


Q_ 

Q0 


T  -T 

bar, initial  bar,  final 


L  bar , initial 


(49) 

(50) 


Time  (sec) 

Q/Qo 

0 

0.000 

1 

0.078 

2 

0.135 

3 

0.197 

4 

0.258 

5 

0.309 

6 

0.355 

8 

0.446 

10 

0.514 

Table  7 :  Q/Q0  ratios  at  given  times. 


7.2  EXPERIMENTAL  VERIFICATION  OF  THE  TRANSIENT  MODEL 

The  transient  model  was  verified  experimentally  by  heating  rail  #  2  (the  optimized  channel 
design)  to  an  initial  unifonn  temperature  and  then  pulsing  water  through  the  channels  for  various 
time  intervals.  In  order  for  water  to  be  pulsed  through  the  rail  for  time  intervals  on  the 
magnitude  of  seconds,  a  bypass  loop  was  routed  into  the  chiller  so  that  the  flow  could  be  diverted 
from  the  bypass  loop  to  the  rail  instead  of  simply  turning  the  flow  through  the  rail  off  and  on  for 
the  specified  time  interval.  In  order  to  divert  the  flow  from  the  bypass  loop  to  the  rail  quickly 
and  for  a  specific  time  interval,  two  solenoid  valves  were  wired  to  a  variable  timing  device.  The 
advantage  of  the  bypass  loop  is  that  it  enables  a  more  unifonn  pressure  and  flow  at  the  beginning 
of  the  time  interval.  Data  was  taken  for  cooling  intervals  from  1.0  to  10.0  seconds. 

The  first  data  set  that  must  be  analyzed  is  the  initial  temperature  profile  along  the  length  of  the 
rail  (Figure  49).  The  temperature  profile  is  not  quite  uniform  due  to  losses  through  the  manifold 
devices,  so  in  order  to  analyze  the  Q/Q0  ratio  without  effects  from  the  heat  loss,  the 
thennocouple  with  the  highest  initial  temperature  was  used  in  the  data  calculations. 


61 


Z,  position  along  length  of  bar  (inches) 


Figure  49:  Temperature  profile  along  rail . 

Once  the  hottest  point  on  the  rail  was  selected  (z  =10),  the  temperature  at  that  point  vs.  time  was 
graphed  (Figure  50).  This  graph  was  used  to  determine  when  the  fast  transient  portion  of  the 
cooling  process  ends  so  that  the  appropriate  final  bar  temperature  could  be  used  in  the  Q/Q0 
calculation. 
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Figure  50:  Temperature  at  thermocouple  10  vs.  time. 


The  experimental  Q/Q0  ratio  is  calculated  as  shown  in  Eq.  50,  except  that  the  final  rail 
temperature  is  determined  from  the  graph  above  rather  than  from  a  weighted  average  of  the 
modeled  temperature  profile.  The  results  show  that  the  bar  is  cooling  down  more  quickly  than 
predicted  (Figure  51).  One  reason  the  Q/Q0  ratio  is  lower  than  predicted  is  that  the  thermal  entry 
region  creates  a  non-uniform  temperature  distribution,  which  causes  heat  to  flow  in  the  z 
direction,  while  the  assumption  used  in  the  model  was  for  conduction  in  only  one  direction. 
Another  reason  that  the  Q/Qa  ratio  is  lower  than  predicted  could  be  that  the  transient  heat  transfer 
coefficient  is  higher  while  the  thermal  boundary  layer  is  developing  within  the  channel. 
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8.0  PROJECT  CONCLUSIONS 

In  order  to  develop  a  long  range  naval  electromagnetic  railgun  that  can  be  fired  in  a  repetitive 
mode,  the  resistive  heating  of  the  rails  in  the  railgun  must  be  considered.  To  determine  the 
optimal  design  for  internal  cooling  channels  in  the  rails,  the  compressive  forces  acting  on  the 
rails  must  also  be  taken  into  account.  In  fact,  these  forces  ultimately  limit  the  optimum  cooling 
channel  design. 

With  the  consideration  of  the  resistive  heating  and  compressive  forces,  an  optimization  routine 
was  created  to  find  the  channel  design  for  optimum  heat  transfer.  This  routine  solves  the  heat 
equation  for  a  fin  with  an  elliptical  boundary  to  find  the  channel  design  that  has  the  lowest 
thennal  resistance  while  meeting  the  defined  structural  constraints. 

The  severity  of  the  stresses  acting  on  the  rails  also  makes  material  selection  an  important  element 
in  the  design.  In  the  cooling  channel  optimization  for  the  full  size  railgun,  a  beryllium  copper 
alloy  was  chosen  because  it  was  the  only  alloy  that  met  the  basic  requirements  for  the  rail 
material  based  on  its  yield  strength  and  thennal  conductivity. 

In  order  to  gain  confidence  in  the  optimization  routine,  the  analytical  expressions  used  to 
represent  the  stress  concentration  factors  and  heat  equation  were  independently  verified.  First 
the  uniaxial  and  biaxial  stress  concentration  factors  used  in  the  optimization  were  confirmed 
using  finite  element  analysis.  For  the  channel  geometries  and  stresses  typical  of  the  railgun 
application,  the  stress  concentration  factors  found  in  the  finite  element  analysis  closely  matched 
the  analytical  stress  concentration  factors.  The  heat  transfer  response  predicted  in  the  cooling 
channel  design  was  verified  through  steady  state  experiments  on  copper  rails  with  three  different 
cooling  channel  designs.  In  this  experiment,  a  heater  was  applied  to  one  surface  of  the  rail,  and 
the  rail  temperature  was  measured  along  the  length  of  the  rail  as  water  was  pumped  through  the 
channels  at  various  flow  rates.  The  temperature  data  collected  was  used  to  calculate  thennal 
resistance,  and  these  values  were  compared  to  the  values  predicted  in  the  optimization  routine. 
The  predicted  values  proved  to  be  in  reasonable  agreement  with  the  experimental  results. 

Because  the  cooling  of  the  rails  in  a  railgun  has  a  transient  nature  in  the  time  between  shots,  the 
transient  thennal  response  of  the  optimized  channel  design  was  also  considered.  An  implicit 
finite  difference  method  was  used  to  predict  the  amount  of  heat  removed  (heat  ratio)  during  a 
given  time  interval,  and  this  quantity  was  then  measured  experimentally  in  order  to  verify  the 
transient  model.  In  the  transient  experiment,  the  rail  was  heated  to  a  uniform  temperature  and 
cooling  water  was  pumped  through  the  channels  for  a  range  of  time  intervals.  From  the 
temperature  data  an  experimental  heat  ratio  was  calculated  for  each  time  interval.  The  predicted 
and  experimental  results  followed  the  same  trend,  but  the  experimental  heat  ratio  was 
consistently  lower  than  what  was  predicted.  This  can  be  attributed  to  the  use  of  one  dimensional 
analysis  and  a  steady  state  heat  transfer  coefficient  in  the  transient  model. 

This  project  provides  an  optimization  routine  with  analytical  expressions  for  heat  transfer  and 
structural  constraints  that  have  been  verified  using  independent  measures.  The  confidence 
afforded  by  the  methodology  employed  in  this  project  indicates  that  this  optimization  routine  is 
suitable  for  further  analysis  of  a  railgun  cooling  channel  design. 
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AREAS  FOR  FURTHER  STUDY 

A  beryllium  copper  alloy  was  used  in  the  full  size  railgun  optimization  out  of  necessity,  but  it  is 
expensive  and  there  are  adverse  health  effects  associated  with  prolonged  exposure  to  beryllium 
compounds  (Haz-Map  2004).  Since  there  is  little  practicality  in  using  such  a  material  in  actual 
railgun  production,  further  research  in  material  selection  is  imperative  for  railgun  design  to 
progress.  Other  potential  alloys  include  chromium-copper  and  GlidCop®  (Appendix  12). 

In  this  study  cooling  channels  were  optimized  for  steady  state  conditions,  so  the  properties 
associated  with  transient  heat  transfer  were  not  considered.  For  example,  thermal  diffusivity,  the 
measure  of  how  fast  heat  is  transferred  through  a  material,  is  not  a  factor  in  the  current  design. 
Since  the  physical  railgun  cooling  problem  has  a  transient  nature,  the  optimization  routine  should 
be  reevaluated  to  consider  transient  performance. 

Many  simplifying  assumptions  were  made  in  this  design — the  most  significant  being  conduction 
in  only  one  dimensional,  steady  state  heating  conditions,  uniform  compressive  loads  and  the 
neglecting  of  any  effects  the  cooling  channels  will  have  on  the  electrical  current  penetration  and 
magnetic  field  distribution.  The  transient  problem  can  be  considered  using  methods  similar  to 
those  used  in  this  study,  but  further  analysis  in  two  or  three  dimensions  introduces  a  level  of 
complexity  that  necessitates  a  finite  element  analysis  approach.  These  computational  studies 
will  need  to  consider  the  coupled  thennal,  mechanical  and  electromagnetic  effects. 
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Appendix  1:  Optimization  Code  (Matlab®) 
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Optimization.m 

%This  is  the  primary  code  of  the  optimization.  Line  13  contains  the  actual 
%function  (f  nincon )  that  does  the  minimization,  fmincon  is  a  function  found  in  the 
%Matlab®  Optimization  Toolbox.  A,  B,  lb  and  ub  are  constraints  on  the  geometric 
%values,  and  pO  is  the  initital  guess.  After  line  13,  the  rest  of  the  code  is  written  to  output 
%values  calculated  with  the  final  solution  (x). 

A=[l  0  0;  0  1  0  ;  0  0  1]; 

B=[0.01  ;  0.012;  0.01]; 
p0=[0.005  ;  0.01  ;  0.005]; 
lb=[0.000 1;0. 000 1;0. 0001]; 
ub=[l;l;l]; 

options  =  optimset(’MaxFunEvals',  500,’MaxIter',  500,’TolCon’,le-50,’TolFun’,le- 
50,’TolX’,le-50); 

[x,fval,exitflag,output]=fmincon(@rfm,pO,  A, B,  [],[], lb, ub,@newstress, options); 

Ncell=3 

x 

mu=x(2)/x(l); 

ta=x(3)/x(l); 

Rtotal=fval/Ncell 

%inputs: 

FS=1; 

yield=l  *3 10*  1000000; 
uniaxial=l  10286000; 

%deH  nitions 

lambda=x(  1 )/ (x(  1  )+x(3)); 
allow=yield/FS; 

%axial  constraint 

nominal=uniaxial  *  (x(  1  )+x(3  ))/x(3) ; 
c  1=  1 . 000+0. 000 *  sqrt(  l/mu)+2. 000*  1/mu; 
c2=-0. 35 1-0.02 1  *sqrt(  l/mu)-2.483*  1/mu; 
c3=  3.621-5.183*sqrt(l/mu)+4.494*l/mu; 
c4=-2.270+5.204*sqrt(  l/mu)-4.0 1 1  *  1/mu; 

Ku=c  1 +c2  *  lambda+c3  *  (lambda)  A2+c4  *  (lambda)A3 

uniaxialmax=Ku*nominal 

%biaxial  constraint 

sigma  1=47124900; 

sigma2=6857 1400; 

width=0.0254; 

effwidth=width-2*x(2); 

sigma  lp=sigma  1  *(width/effwidth) 

sigma2p=sigma2*((x(l)+x(3))/x(3)) 

biaxialstressratio=sigma  1  p/sigma2p 

if  sigma  lp/sigma2p  <  1/mu; 
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Kb=-sigmalp/sigma2p+  1+2/mu 
biaxialmax=sigma2p  *  ( 1  +2/mu)-sigma  1  p 
’middle  L' 

else  Kb=(sigma  lp/sigma2p)*(  l+2*mu)- 1 
biaxialmax=(sigma  1  p  *  ( 1  +2  *mu)-sigma2p) 

’positive  SL' 
end 

%Reynolds  number  for  solution: 

rho=997; 

mu=855e-6; 

Aellipse=pi*x(l)*x(2); 

P=pi*(3  *(x(  1  )+x(2))-sqrt((3  *x(  1  )+x(2))*(x(  1  )+3*x(2)))); 

Dh=4*Aellipse/P 

Um=feval(@Umm,x) 

Re=rho*Dh*Um/mu 

%heat  transfer  coefficient  for  solution: 

h=feval(@hxcoef,x) 

exitflag 

rfin.m 

%This  function  defines  the  fin  resistance  (the  quantity  being  minimized),  which  involves 
solving  the 

%boundary  value  problem  or  the  energy  equation  for  a  variable  area  fin. 

%The  function  bvp4c  is  used  to  solve  this  equation  (testg)  with  its  boundary 
%conditions  (testbcg).  The  definition  of  rfin  is  the  last  line  of  this 
%code  (line  13). 

function  x=rfin(g); 
k=388; 

sohnit=bvpinit(linspace(-0.99*g(2),0.99*g(2),10),[l,0]); 

sol=bvp4c(@testg,@testbcg,solinit,[],g); 

y=linspace(-0.99*g(2),0.99*g(2)); 

T=deval(y,sol); 
x=  l/(-k*(g(  1  )+g(3))*T(2)); 

%hold  on 
%plot(y,T(2,:)) 


testg.m 

%This  function  contains  the  energy  equation  for  a  fin  with  elliptical 
%boundaries.  It  is  defined  in  terms  of  y,  T  and  g,  where  g  is  an  array  if 
%(a,b,t).  The  heat  transder  coefficient  is  also  dependent  on  g,  and  it  is 
%defined  here  as  a  function  handle  that  calls  the  function  hxcoef.m 
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function  dTdy=testg(y,T,g); 

h=feval(@hxcoef,g); 

k=388; 

xstar=g(  1  )+g(3)-sqrt((  1  -yA2/ g(2)A2)*g(  1  )A2); 

dTdy=[T  (2)  ;(-(g(  1 )/ g(2))  *y  *  T(2)/ sqrt(g(2)A2-yA2)+(h/k)*  T  ( 1 )  *  sqrt(  1 +(- 
y  *(g(  1 V  g(2))/ sqrt(g(2)A2-yA2))A2))/xstar] ; 

testbcb.m 

%This  function  contains  the  boundary  condidtions  for  the  function  in 
%testg.m 

function  res=testbc(Ta,Tb,g) 
res=[Ta(  1 )- 1  ;Tb(2)] ; 


hxcoef.m 

%The  heat  transfer  coefficient  is  a  function  of  many  variables,  many  of  which  are 
%affected  by  the  geometry  of  the  channel  (g).  The  fluid  properties  are  assumed  to 
%be  constant  and  taken  at  room  temperature.  The  mean  fluid  velocity  (Um)  must  be 
%detennined  by  a  separate  function  (Umm.m). 

function  x=hxcoef(g); 

Um=feval(@Umm,g); 

rho=997; 

mu=0. 0008905; 

Cp=4183; 

kf=0.5948; 

Aellipse=pi*g(l)*g(2); 

P=pi*(3  *(g(  1  )+g(2))-sqrt((3  *g(  1  )+g(2))*(g(  1  )+3*g(2)))); 

Dh=4*Aellipse/P; 

Re=rho*Dh*Um/mu; 

Pr=Cp*mu/kf; 
f=(0.79*log(Re)- 1 .64)A(-2); 

Nu=((f/8)*(Re-1000)*Pr)/(l+12.7*sqrt(f/8)*(PrA(2/3)-l)); 

x=Nu*kf/Dh; 

pumpconstraint.m 

%The  pump  constraint  in  this  case  is  ?  Watts/meter. 

%pumppower=10/Ncell  (per  unit  length) 

%pumppower=pressuredrop*Aellipse*Um  (per  unit  length) 
%10/Ncell=pressuredrop*Aellipse*Um  (per  unit  length) 

%in  order  to  find  the  value  of  Um,  the  equation  is  manipulated  to  so  that 
%the  zero  of  the  equation  will  be  the  value  of  Um: 
%x=pressuredrop*Aellipse*Um-?/Ncell 


function  x=pumpconstraint(Um,g); 

Ncell=3; 

rho=997; 

mu=855e-6; 

Aellipse=pi*g(l)*g(2); 

P=pi*  (3  *  (g(  1  )+g(2))-sqrt((3  *  g(  1  )+g(2))  *(g(  1  )+3  *  g(2))));  %Perimeter% 
Dh=4  *  Aellipse/P; 

Re=rho*Dh*Um/mu; 

f=(0.79 1  *log(abs(Re))-1.64)A(-2); 

prcssurcdrop=P  rho  *Um  A2/ (2  *  Dh) ; 

x=pressuredrop*  Aellipse*Um- 1 .9 1 84 1 50576964/Ncell; 

Umm.m 

%This  mean  fluid  velocity  is  determined  by  finding  the  zero  of  the  pump 
%constraint  equation  (pumpconstraint.m) 

function  x=Umm(g); 
x=biseet(’pumpconstraint’,0.0 1 ,30,g); 

bisect.m 

function  root  =  bisect  (  f,  a,  w,  g) 

fa  =  feval  (  f,  a,  g); 
fw  =  feval  (  f,  w,  g); 

if  (  sign  (  fa  )  ==  sign  (  fw  )  ) 

’BISECT  -  Fatal  error!’ 

’  [A,w]  is  not  a  change-of-sign  interval.’ 
root  =  0; 
return 
end 

format  long 

while  (  abs  (  w  -  a  )  >  0.000001  ) 

c  =  (a  +  w)/2; 
fc  =  feval  (  f,  c,  g  ); 

[  a,  c,  w; 
fa,  fc,  fw  ]; 

if  (  fc  ==  0  ) 
root  =  c; 
return 

elseif  (  sign  (  fc  )  ==  sign  (  fa  )  ) 
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a  =  c; 
fa  =  fc; 

elseif  (  sign  (  fc  )  ==  sign  (  fw  )  ) 
w  =  c; 
fw  =  fc; 
end 

end 

root  =  0.5  *  (  a  +  w  ); 


newstress.m 

%This  file  contains  all  of  the  structural  constraints  that  are  place  on  the 
%optimization. 

function  [c,ceq]=newstress(g); 

%inputs: 

FS=1; 

yield=l  *3 10*  1000000; 
allow=yield/FS; 

Ncell=3; 

height=0.0508; 

cellwidth=height/Ncell; 

uniaxial=l  10286000; 

sigma  1=47124900; 

sigma2=6857 1400; 

width=0.0254; 

effwidth=width-2  *  g(2) ; 

sigma  lp=sigma  1  *(width/effwidth); 

sigma2p=sigma2*((g(l)+g(3))/g(3)); 

%definitions 
mu=g(2)/g(l); 
lambda=g(  1  )/(g(  1  )+g(3)); 

%axial  constraint 

nominal=uniaxial*(g(l)+g(3))/g(3); 
cl=  1. 000+0. 000*sqrt(l/mu)+2. 000*  1/mu; 
c2=-0. 35 1-0.02 1  *sqrt(  l/mu)-2.483*  1/mu; 
c3=  3.621-5.183*sqrt(l/mu)+4.494*l/mu; 
c4=-2. 270+5 .204*sqrt(  l/mu)-4.0 1 1  *  1/mu; 

Kc=c  1 +c2  *  lambda+c3  *  (lambda)A2+c4  *  (lambda)  A3 ; 

%biaxial  constraint 
if  sigma  lp/sigma2p  <  1/mu; 

K=-sigma  1  p  /  sigma2p+ 1  +2/ mu; 
sigmamax=sigma2p  *  ( 1  +2/mu)-sigma  1  p ; 

’middle  line  used'; 


else  K=(sigmalp/sigma2p)*(l+2*mu)-l ; 
sigmamax=(sigma  1  p  *  ( 1 +2*mu)-sigma2p); 
’positive  slope  line  used'; 
end 

%minimum  effective  height  constraint 

sigmaarea=uniaxialHsheight/(Ncell*2*g(3)); 

c=[Kc*nominal-allow;sigmamax-allow;0]; 

%mineffh-2*g(3)*Ncell 

%unit  cell  constraint 

ceq=[g(l)+g(3)-cellwidth/2;0;0]; 
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Appendix  3:  Average  Surface  Temperature  Calculations  (Matlab®) 


Conclusioncalcs.m 

%  final  calculations  for  conclusion 

for  i=  1:1:50 
time(i)=i*.5; 
end 

Heat=15000000;  %Joules 
width=0.135;  %  meters 
length=10;  %  meters 
Area=length*width; 
q_flux=Heat./(time*Area);  %Watts 
d=0.020;  %meters 
k=215;  %W/m-K 

R_doubleprime_conduction=d/k; 

R_cell_per_Length=0.0007 18707497568427;  %for  8  channels  in  the  full  size 
R_doubleprime_fin=R_cell_per_Length*width; 

R _ doubleprime_total=R_doubleprime_fin+R_doubleprime_conduction; 

T_inf=293; 

Delta_T=q_flux*R _ doubleprime  total;  %Kelvin  or  C 

T_surface=T_inf+Delta_T-273 ;  %Celsius 

hold  on 

plot(time,T_surface, ’-black') 
plot(time(  1 0),T_surface(  1 0),'oblack') 
plot(time(20),T_surface(20),'oblack') 
axis([0  25  0  4500]) 
xlabel(’time  (seconds)') 
ylabel(’Temperature  (C)') 


Appendix  4:  Finite  Element  Analysis:  Uniaxial  Stress 
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Notes: 

— Maximum/Minimum  Principal  Stress  is  the  Stress  Concentration  Factor  from  l-DEAS  because  the 
model  was  loaded  with  a  force  of  1 .0 

— Analytical  stress  concentration  factor  comes  from  Young  and  Budynas  (2002). 

—  %  error  refers  to  the  diffference  between  the  l-DEAS  and  analytical  stress  concentration  factors 
— See  page  5  of  this  appendix  for  plate  dimensions  and  finite  element  model  parameters 

1 .  Verify  that  a  circular  hole  in  a  plate  loaded  in  tension  and  compression  have  the  same  stress 

concentration  factors 

Analytical  Stress 


l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Concentration  Factor: 

quarter  piece  tension 

tension 

tension  soln  setl 

4.344 

Results: 

max  stress  in  x 

4.349 

max  stress  in  y 

0.588 

max  stress  in  z 

1.262 

%  error 

Maximum  Principal  stress 

4.349 

-0.12 

Von  Mises  stress 

3.874 

l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Analytical  Stress 
Concentration  Factor: 

quarter  piece  compression 

compression 

compression  soln  setl 

4.344 

Results: 


max  stress  in  x 

-4.349 

max  stress  in  y 

-0.588 

max  stress  in  z 

-1.262 

Minimum  Principal  stress 

-4.349 

Von  Mises  stress 

3.874 

Conclusion:  Tension  and  Compression  have  stress  the  same  stress  concentration  factors 
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2.  Verify  mesh  independence  of  solution 


Analytical  Stress 


l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Concentration  Factor: 

quarter  piece  compression 

compression  2 

compression  soln  set2 

4.344 

Results: 

max  stress  in  x 

-4.355 

max  stress  in  y 

-0.588 

max  stress  in  z 

-1.26 

%  error 

Minimum  Principal  stress 

-4.355 

-0.26 

Von  Mises  stress 

3.88 

l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Analytical  Stress 
Concentration  Factor: 

quarter  piece  compression 

compression  2 

compression  soln  set2 

4.344 

Results: 

(2nd  time)* 

max  stress  in  x 

-4.3458 

%  error 

Minimum  Principal  stress 

-4.3467 

-0.06 

*  The  number  of  elements  was  increased  in  order  to  verify  mesh  independence 
Conclusion:  Mesh  independence  achieved 


3.  Verify  element  shape  independence  in  solution 

Analytical  Stress 

l-DEAS  Part  name:  FE  model  name:  l-DEAS  solution  name:  Concentration  Factor: 


quarter  piece  compression  [compression  3  [compression  soln  set3 


4.344 


Results: 


max  stress  in  x 

-4.333 

%  error 

Minimum  Principal  stress 

-4.334 

0.24 

l-DEAS  Part  name: 

FE  model  name: 

Analytical  Stress 

l-DEAS  solution  name:  Concentration  Factor: 

quarter  piece  compression 

compression  4 

compression  soln  set4  4.344 

Results: 

max  stress  in  x 

-4.3557 

%  error 

Minimum  Principal  stress 

-4.3569 

-0.30 

Conclusion:  the  solution  is  practically  independent  of  element  shape,  but  the  three  point  triangle 
element  (FE  Model:  compression  3)  gave  slightly  better  results 
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4.  Verify  that  the  numerical  solution  matches  the  analytical  solution  for  a  range  of  plate  sizes 

(hole  is  the  same  size  and  relative  location) 


I-DEAS  Part  name: 

FE  model  name: 

Analytical  Stress 

l-DEAS  solution  name:  Concentration  Factor: 

quarter  2 

compression  5 

compression  soln  set5  4.344 

Results: 

max  stress  in  x 

-4.597 

%  error 

Minimum  Principal  stress 

-4.598 

-5.84 

l-DEAS  Part  name: 

FE  model  name: 

Analytical  Stress 

l-DEAS  solution  name:  Concentration  Factor: 

quarter  piece  compression 

compression  3.5 

compression  soln  set3.5  4.344 

Results: 

max  stress  in  x 

-4.3228 

%  error 

Minimum  Principal  stress 

-4.3228 

0.49 

Note:  these  cases  were  run  to  investigate  the  effect  of  overall  plate  length 

Conclusion:  the  length  of  the  plate  in  the  direction  of  the  force  does  affect  the  stress  concentration 
factor,  even  though  the  equation  from  the  literature  does  not  account  for  this  effect. 
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5.  Verify  that  the  numerical  solution  matches  the  analytical  solution  for  an  elliptical  hole 

Analytical  Stress 


I-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Concentration  Factor: 

quarter  3 

compression  6 

compression  soln  set  6 

2.302 

Results: 

max  stress  in  x 

-2.260 

%  error 

Minimum  Principal  stress 

-2.261 

1.77 

l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Analytical  Stress 
Concentration  Factor: 

quarter  4 

compression  7 

compression  soln  set7 

1.99 

Results: 

max  stress  in  x 

-2.01 

%  error 

Minimum  Principal  stress 

-2.01 

-1.01 

l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Analytical  Stress 
Concentration  Factor: 

quarter  5 

compression  8 

compression  soln  set8 

6.565 

Results: 

max  stress  in  x 

-6.03 

%  error 

Minimum  Principal  stress 

-6.04 

8.00 

l-DEAS  Part  name: 

FE  model  name: 

l-DEAS  solution  name: 

Analytical  Stress 
Concentration  Factor: 

quarter  6 

compression  9 

compression  soln  set9 

5.545 

Results: 

max  stress  in  x 

-5.27 

%  error 

Minimum  Principal  stress 

-5.29 

4.60 

Conclusion: 

For  an  ellipse  the  correlation  is  better  when  the  aspect  ratio  (beta)  is  greater  than  one. 
However,  for  an  aspect  ratio  less  than  1 ,  the  analytical  solution  is  conservative. 


Part  Dimensions 
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Notes: 

—  -i  is  acting  in  the  y  direction  (up  and  down),  2  acts  in  the  x  direction 

—  a  is  in  the  y  direction,  b  is  in  the  x  direction 
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AMERj  CAN_  ML  CR  O Industries,  Inc. 

440  C  Ramsey  Avenue,  Chambersburg,  PA  17201 
Tel:  (800)  558-2058  Fax:  (717)  261-9161  e-mail:  sales@americanmicroinc.com 


Material  Specifications  for  G10/FR4 


Properties 

NEMA  grade  reinforcement-resin  binder,  glass-epoxies 

Tensile  Strength 

Lengthwise,  PSI 

40,000 

Crosswise,  PSI 

35,000 

Compressive  Strength 

Flatwise,  PSI 

60,000 

Edgewise,  PSI 

35,000 

Flexural  Strength 

Lengthwise,  PSI 

55,000 

Edgewise,  PSI 

45,000 

Modulus  of  Elasticity  in  Flex 

Lengthwise,  PSI  x  10  6 

2.7 

Crosswise,  PSI  x  10  6 

2.2 

Shear  Strength,  PSI 

19,000 

IZOD  Impact 

Flatwise,  ft  lb  per  inch  of  notch 

7 

Edgewise,  ft  lb  per  inch  of  notch 

5.5 

Rockwell  Hardness  M  Scale 

110 

Specific  Gravity 

1.82 

Coefficient  of  Thermal  Expansion 

cm/cm  deg  C  x  10  5 

.9 

Water  Absorption 

0.62”  thick,  %  per  24  hrs 

0.25 

0.125”  thick,  %  per  24  hrs 

0.15 

0.500”  tick,  %  per  24  hrs 

0.10 

Dielectric  Strength,  volt/mil 

0.062”  thick  (perpendicular  to  laminations;  short) 

500 

0.125”  thick  (perpendicular  to  laminations;  short) 

400 

Dissipation  Factor 

condition  A,  1  megacycle 

0.025 

Dielectric  Constant 

condition  A,  1  megacycle 

5.2 

Insulation  Resistance 

Condition:  96  hrs  at  90%  relative  humidity  (in  mega  ohms) 

200,000 

Flame  Resistance 

Underwriter  Labs,  Classification 

94V-0 

Bond  Strength 

in  lbs 

2,000 

Max  Continuous  Operating  Temperature 

Approximate  degrees  °F 

285 

The  values  presented  are  typical  and  not  intended  for  specification  purposes.  This  information  is  provided  without  warranty,  representation,  inducement  or  license  of  any 
kind,  including,  but  not  limited  to  implied  warranties  of  merchantability  and  fitness  for  a  particular  use  or  purpose,  except  that  it  is  accurate  to  the  best  of  American  Micro 
Industries'  knowledge  or  obtained  from  sources  believed  by  American  Micro  Industries  to  be  accurate,  and  American  Micro  Industries  does  not  assume  any  legal 
responsibility  for  the  use  or  reliance  upon  same,  as  well  as  any  typographical  errors.  Users  of  our  product(s)  are  encouraged  to  conduct  their  own  tests  for  suitability  and 
conformance. 


Rev.  5/17/04 


Appendix  8:  Thermocouple  Calibration 


icebath2.m 

%This  file  contains  the  data  from  the  thermocouple  calibration,  and  running  this  produces  the 
%matrices  'ice'  and  'boil'  which  are  used  in  the  steady  state  calculations.  Consequently,  this  file 
%must  be  run  before  the  files  containing  the  steady  state  %calculations  must  be  run 


x=[0;l;2;3;4;5;6;7;8;9;10;ll;12;13;14;15;16;17;18]; 


ice=[0;0.217595 194805 19487 ; 
0.12027805194805201; 
0.035655194805194815; 
0.1989222077922078; 
;0.332133 116883 11703; 


0.14695818181818177; 

0.2511764935064935; 

0.3833594805194806; 

0.02924597402597404; 


0.2727576623376624; 

0.09331064935064941; 

0.04666545454545456; 

0.10564740259740257 


0.24440649350649354;  0.45613961038961043; 


0.3643548051948051; 


0.42522597402597395  ;0]; 


figure(2) 

plot(x,ice,'oblue') 

title('Thermocouple  reading  in  an  icebath  2  vs.  thermocouple  number')  % 
xlabel  ('thermocouple  number') 
ylabel  ('Temperature  (C)')  % 


boil=[100;  100.0648827777778;  100.066301 11111112;  100.30158266666662; 

100.2308548888889;  100.19585133333337;  100.04096233333333; 

100.00681766666662;  100.38428022222222;  100.06214277777781; 

100.20797055555555;  100.13512822222219;  100.24105055555557; 

100.2351358888889;  100.36998788888886;  100.45128022222222; 

100.32367366666666;  100.4432161 1 1 1 1 106;  100]; 


figure(3) 

plot(x,boil,'ored') 

title('Thermocouple  reading  in  an  boiling  distilled  water  vs.  thermocouple  number')  % 
xlabel  ('thermocouple  number') 
ylabel  ('Temperature  (C)')  % 

difference=boil-ice 


figure(4) 

plot(x, difference, 'omagenta') 

title('Temperature  difference  (boil- ice)  vs.  thermocouple  number')  % 
xlabel  ('thermocouple  number') 
ylabel  ('Temperature  (C)')  % 

%Note:  points  0  and  18  were  not  actually  measured,  they  are  place  holders 
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This  is  a  sample  code  used  to  manipulate  the  steady  state  data.  The  subsequent  codes  for 
various  flow  rates  and  channel  geometries  are  identical  except  for  the  rail  temperatures 
(LPM2)  and  water  temperatures  (waterIN2LPM  and  waterOUT2LPM). 

iLPM2data.m 

%in  order  to  get  the  legend  to  work  correctly,  all  the  LPM2data,  etc  files 
%need  to  be  run  in  order  by  increasing  flow  rate, 
hold  on 

axis([0  18  0  12]) 

title(’Temperature  Difference  vs.  Length  Along  Bar')  % 
xlabel  (’X  (in)') 

ylabel  (’Temperature  Difference  (C)')  % 
x=[0;l;2;3;4;5;6;7;8;9;10;ll;12;13;14;15;16;17;18]; 

LPM2=[0;23.99;25.527;27.085;27.89;28.886;29.219;29.409;30.168;29.737;30.218;29.87 

;30.176;30.582;0;30.47;30.565;29.774;0;]; 

waterIN2LPM=[19.640151 1 1;19.53934189;19. 50703789]; 
waterOUT2LPM=[22.0249581 1;22. 02905856;22. 01 340467]; 

IN=sum(waterIN2LPM)/3 ; 

OUT=sum(  waterOUT  2LPM)/3 ; 
deltay=(OUT -IN) ; 
deltax=18; 
m=deltay/deltax; 

y=m*x; 

plot(x,y, 'black') 

theta2=(LPM2-ice)  *100./ (boil-ice)-IN ; 
plot(x,theta2,'oblack') 

length=  18 *0.0254;  %meters 

qPerLength=(1000/6)/length;  %Watts/meter%  there  are  6  half-unit  cells  in  the  bar 
for  i=  1:19 

deltaT(i)=theta2(i)-y(i); 

uaPerLength(i)=qPerLength/deltaT(i); 

rfinPerLength(i)=deltaT(i)/qPerLength; 

end 

UAPerLength=uaPerLength' 

R_finPerLength=rfinPerLength’ 


%'Note:  units  for  UA  are  W/m-K,  units  for  R_fin  are  K-m/W' 


%'the  negative  values  do  not  correlate  with  data,  they  are  only  being  used  as 
placeholders.' 
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R_sum=R_finPerLength(8)+R_finPerLength(9)+R_finPerLength(  1 0)+R_finPerLength(  1 
1  )+R_finPerLength(  1 2)+R_finPerLength(  1 3)+R_finPerLcngth(  1 4); 
R_fin_PerLength_AVERAGE=R_sum/7 


Appendix  10:  Steady  State  Data  Calculations 
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Appendix  10:  Steady  State  Data  Calculations 
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Appendix  11:  Transient  Model 


middletemp.m 

%Inputs  from  experiment: 

%  time  h  —  based  on  flow  rate  and  channel  geometry 

clear  all 

for  m=l:l:10 

Time=[l;2;3;4;5;6;7;8;9;10];%  The  time  is  entered  this  way  instead  of  just  through 
for  loop  so  I  can  put  in  times  that  are  not  nice  clean  numbers  like  here 
time(m)=Time(m); 

a=0. 00328086088597;  %  from  optimization  with  experimental  parameters 
b=0. 00605032567597;  %  from  optimization  with  experimental  parameters 
t=0. 005 1858057807;  %  from  optimization  with  experimental  parameters 
cell_dim_y=a+t; 

cell_dim_x=  1*0.0254;  %1  inch  (length  of  experimental  rail  cross  section) 
d=0.5Hs(cell_dim_x-2*b);  %  distance  between  edge  of  rail  and  channel  (x  direction) 
rho_bar=8936;  %kg/m,  from  EES,  T=T_inf 
Cp_bar=383;  %J/kg-K,  from  EES,  T=T_inf 
k_bar=388;  %W/m-K,  from  EES,  T=T_inf 

h_water=7.249582579872966e+003;  %%%%%  flow  rate  affects  h — from  to  hxcoeff. 
b_node_quantity=  1 00; 
dx_b=2  *  b  /  (b_node_quantity- 1 ) ; 
dt=0.0005; 

%  time  =10;  %%%  use  this  instead  of  initial  for  loop  (which  creates  a 
%  matrix  of  time  values)  when  you  want  a  solution  for  only  one  time 
%  interval 

number_of_timesteps=round(time(m)/dt); 

T_inf=293;  %%%%Kelvin 

Fo_b=(k_bar*dt)/(dx_bA2*rho_bar*Cp_bar);  %  Fourier  number 

d_node_quantity=floor(d/dx_b); 

%floor  function  rounds  down  to  the  nearest  integer 
extra_node_length=d-d_node_quantity  *  dx_b ; 

%this  accounts  for  the  extra  length  that  will  be  added  to  the  elements 
%on  each  end  of  the  unit  cell.  This  extra  length  is  a  product  of 
%placing  a  node  at  the  endge  of  the  ellipse  and  maintaining  a  uniform 
%dx  across  the  unit  cell 

total_node_quantity=2  *  dnodequantity+bnodequantity ; 
dx=dx_b; 

for  i=l :  1  :total_node_quantity 
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xj(l)  =  0; 

xj‘(i)  =  (dx*(i-l))+extra_node_length; 
xj'(total_node_quantity)  =  (dx*(i-l))+2*extra_node_length; 

%xj  and  yj  values  are  midway  between  nodes,  to  the  right  of  the  node 
end 

for  i=l :  1  :total_node_quantity-l 

%this  for  loop  defines  the  x  and  y  coordinates  of  the  ellipse  across 
%the  whole  unit  cell 

if  i  <=  dnodequantity 
x(i)=0; 

y_ellipse(i)  =0; 

elseif  i  >=  (dnodequantity+bnodequantity) 
x(i)=0; 

y_ellipse(i)  =  0; 


else 

x(i)  =  (xj(i)+xj(i+l))/2  -  (b+d)  ; 

y_ellipse(i)  =  sqrt(aA2*(l-x(i)A2/bA2));  %to  check,  plot(xj,y_ellipse,'o')  axis([0  0.06  - 
0.06  0.06]) 

end 

yj  =  cell_dim_y  -  y_ellipse; 
end 

yj(totalnodequantity)  =  celldimy; 
for  i=l :  1  :total_node_quantity 

%  this  for  loop  defines  the  length  of  the  channel  boundary  for  each 
%  element 

if  i  <  d  node  quantity 
LinS(i)=0; 

elseif  i  ==  d_node_quantity+l 

LinS(i)  =  sqrt  ((x(i)— b)A2  +  (y_ellipse(i)-0)A2); 

elseif  i  ==  dnodequantity+bnodequantity 

LinS(i)  =  sqrt  ((b-x(i-l))A2  +  (y_ellipse(i- 1  )-0)A2); 

elseif  i  >  d  node  quantity+b  node  quantity 
LinS(i)=0; 


else 


LinS(i)  =  sqrt  ((x(i)-x(i-l))A2  +  (y_ellipse(i)-y_ellipse(i-l))A2); 

%note:  for  verification,  we  checked  that  sum(LinS)=Perimeter  of  the 
%ellipse/2 
end 
end 

yf(l)  =  cell_dim_y; 
yf(total_node_quantity)  =  celldimy; 

for  i=2: 1  :total_node_quantity- 1 ; 

%this  for  loop  defines  the  values  of  yf,  the  average  height  of  an  element 
yf(i)=(yj(i)+yj(i-i))/2; 
end 

yf(d_node_quantity+l)  =  3/4*yj(d_node_quantity)  +  l/4*yj‘(d_node_quantity+l); 
yf(d_node_quantity+b_node_quantity)  =  l/4*yj(d_node_quantity+b_node_quantity- 1) 
3/4  *  yj(d_node_quantity+b_node_quantity); 

%these  yf  values  are  for  the  elements  whose  node  falls  at  tne  boundary  of 
%the  channel 

for  i=l :  1  :total_node_quantity 

T(i)=273+300;  %%%%  This  is  the  intial  temperature  of  the  unit  cell 
end 

T_initial=T(i); 

for  n=l :  1  :number_of_timesteps; 

%for  each  time  step,  the  matrix  will  be  generated  and  the 
%temperatures  across  the  unit  cell  will  be  solved  for 
for  i=l :  1  :total_node_quantity; 

%this  for  loop  defines  column  vectors  of  the  coefficients  a,  b,  c 
%and  d  which  are  tenns  in  the  heat  equation 
if  i  ==  1 
a(i)  =  0; 
else 

a(i)  =  -F  o_b  *yj  (i- 1  )/((yj  (i)+yj  (i- 1  ))/2); 

%  a(i)  represents  the  energy  coming  from  the  left 
%  hand  side 
end 

b(i)  =  1  +  2*Fo_b  +  h_water*LinS(i)*dt/(rho_barH5Cp_bar*dx*(yf(i))); 

%  b(i)  represents  the  energy  coming  from  the  node  itself 
if  i  ==  totalnodequantity 
c(i)  =  0; 
elseif  i==  1 

c(l)  =  -Fo_b*yj(l)/yf(l); 
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else 

c(i)  =  -F  o_b  *yj  (i)/ ((yj  (i)+yj  (i- 1  ))/2); 

%  c(i)  represents  the  energy  coming  from  the  right  hand  side 
end 

d(i)  =  T(i)  +  h_water*LinS(i)*dt*T_inf/(rho_bar*Cp_bar*dx*(yf(i))); 

%  d(i)  represents  the  energy  at  timestep  n  (as  opposed  to  the 
%  other  column  vectors  which  represent  the  energy  at  timestep 
%  n+1) 
end 

a(totalnodequantity)  =  -Fo_b*yj(total_nodc_quantity- 1  )/yf(total_node_quantity); 
b(totalnodequantity)  =  1  +  Fo_b*(yj(total_node_quantity- 
l)/yf(total_node_quantity))  + 

h_water*LinS(total_node_quantity)*dt/(rho_bar*Cp_bar*dx*yf(total_node_quantity)); 
%this  changes  the  right  (back)  edge  node  equation,  we  say  that  the 
%right  edge  is  insulated, 

c(l)  =  -Fo_b*yj(l)/yf(l); 

b(l)=  1  +  Fo_b*(yj'(l)/yf(l))  +  hjwater*LinS(l)Mt/(rho_bar*Cp_bar*dx*yf(l)); 
%this  changes  the  left  (front)  edge  node  equation,  we  say  that  the  left 
%edge  is  insulated, 

T  =  Trimatrix(total_node_quantity,a,b,c,d); 

%Trimatrix  is  a  sparse  matrix  solver 

!’• 

1  5 

for  i=l :  1  :total_node_quantity 
T_final=T ; 

Syf=sum(yf); 

yfT(i)=yf(i)*T(i); 

Syfr=sum(yfT); 

T_final_average=SyfT/Syf; 

%This  is  a  final  temperature  computed  as  a  weighted  average  based 
%on  the  volume  (which  reduced  to  yf)  of  each  node 
end 

if  n==number_of_timesteps 

Q_ratio(m)=(T_initial-T_final_average)/(T_initial-T_inf); 

%  This  calculates  the  dimensionless  heat  ratio,  Q/Qo,  the  ratio  of 
%  the  heat  removed  to  the  total  amount  of  heat  that  could  be  removed 
end 

hold  all 

k=l :  1  :total_node_quantity; 
figure(m) 


if  n==round(number_of_timesteps/ 120) 
plot(k,T,’red') 
end 

if  n==round(number_of_timesteps/20) 
plot(k,T, ’magenta’) 
end 

if  n==round(number_of_timesteps/ 6) 
plot(k,T,’yellow’) 
end 

if  n==round(number_of_timesteps/3) 
plot(k,T,’green’) 
end 

if  n==round(number_of_timesteps/l  .5) 
plot(k,T,’cyan’) 
end 

if  n==number_of_timesteps 
plot(k,T,’blue’) 
end 

end 

axis([0  total_node_quantity  T  inf  573]) 

%title(’Temperature  vs.  Node  Number’, ’fontsize’, 20) 
xlabel(’Node  Number  (x-direction)’) 
ylabel(’Temperature  (Kelvin)’) 

%N.B.:  This  legend  is  only  accuurate  for  the  first  timestep 
%legend(’25  timesteps',’150  timesteps','500  timesteps',  T000  timesteps', 
%'2000  timesteps',  ’3000  timesteps', 'location', 'EastOutside') 
end 

Qratio' 

%  specific  legends  for  figures  1  and  10: 

%  figure(l)  legend(’0.0083  seconds', '0.040  seconds',’0.167  seconds', '0.333 
%  seconds’, '0.667  seconds', T. 000  seconds', 'location', 'Eastoutside') 

%  ylabel(’Temperature  (Kelvin)')  xlabel(’Node  Number  (x-direction)') 

%  titleO  ’) 

% 

%  figure)  10)  legend(’0.083  seconds', '0.40  seconds', T. 67  seconds', '3. 33 
%  seconds', '6. 67  seconds’, TO. 0  seconds’, ’location’, 'Eastoutside') 

%  ylabel(’Temperature  (Kelvin)')  xlabel(’Node  Number  (x-direction)') 

%  title)'  ’) 


Trimatrix.m 

function  T  =  Trhnatrix(n,A,B,C,RHS) 
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E(n)  =  0; 
F(n)  =  0; 
T(n)  =  0; 


for  i=l:l:n 

E(i  +  1)  =  -(1  /  ((B(i)  +  A(i)  *  E(i))))  *  C(i); 

F(i  +  1)  =  (1  /  ((B(i)  +  A(i)  *  E(i))))  *  (RHS(i)  -  A(i)  *  F(i)); 
end 

% —  Solution  for  the  outer  boundary - 

T(n)  =  F(n  +  1); 

%—  Solution  of  the  remaining  by  back  substitution - 

for  i  =  (n- 1):- 1 : 1 

T(i)  =  E(i  +  1)  *  T(i  +  1)  +  F(i  +  1); 
end 


Appendix  12:  Material  Properties  Search:  Possible  Rail  Materials 
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Reference:  Matweb.com 

Materials  listed  in  order  of  decreasing  electrical  resistivity 


